Mercurial > hg > octave-lyh
annotate liboctave/dbleLU.cc @ 11188:4cb1522e4d0f
Use function handle as input to cellfun,
rather than quoted function name or anonymous function wrapper.
| author | Rik <octave@nomad.inbox5.com> |
|---|---|
| date | Wed, 03 Nov 2010 17:20:56 -0700 |
| parents | 12884915a8e4 |
| children | 141b3fb5cef7 |
| rev | line source |
|---|---|
| 457 | 1 /* |
| 2 | |
| 7017 | 3 Copyright (C) 1994, 1995, 1996, 1997, 1999, 2002, 2003, 2004, 2005, |
| 8920 | 4 2007, 2008 John W. Eaton |
| 9708 | 5 Copyright (C) 2009 VZLU Prague, a.s. |
| 457 | 6 |
| 7 This file is part of Octave. | |
| 8 | |
| 9 Octave is free software; you can redistribute it and/or modify it | |
| 10 under the terms of the GNU General Public License as published by the | |
| 7016 | 11 Free Software Foundation; either version 3 of the License, or (at your |
| 12 option) any later version. | |
| 457 | 13 |
| 14 Octave is distributed in the hope that it will be useful, but WITHOUT | |
| 15 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or | |
| 16 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License | |
| 17 for more details. | |
| 18 | |
| 19 You should have received a copy of the GNU General Public License | |
| 7016 | 20 along with Octave; see the file COPYING. If not, see |
| 21 <http://www.gnu.org/licenses/>. | |
| 457 | 22 |
| 23 */ | |
| 24 | |
| 25 #ifdef HAVE_CONFIG_H | |
| 1192 | 26 #include <config.h> |
| 457 | 27 #endif |
| 28 | |
| 29 #include "dbleLU.h" | |
| 1847 | 30 #include "f77-fcn.h" |
| 457 | 31 #include "lo-error.h" |
| 9708 | 32 #include "oct-locbuf.h" |
| 33 #include "dColVector.h" | |
| 457 | 34 |
| 1992 | 35 // Instantiate the base LU class for the types we need. |
| 36 | |
| 37 #include <base-lu.h> | |
| 38 #include <base-lu.cc> | |
| 39 | |
|
8367
445d27d79f4e
support permutation matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7482
diff
changeset
|
40 template class base_lu <Matrix>; |
| 1992 | 41 |
| 42 // Define the constructor for this particular derivation. | |
| 43 | |
| 457 | 44 extern "C" |
| 45 { | |
| 4552 | 46 F77_RET_T |
| 5275 | 47 F77_FUNC (dgetrf, DGETRF) (const octave_idx_type&, const octave_idx_type&, double*, |
|
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
48 const octave_idx_type&, octave_idx_type*, octave_idx_type&); |
| 9708 | 49 |
| 50 #ifdef HAVE_QRUPDATE_LUU | |
| 51 F77_RET_T | |
| 52 F77_FUNC (dlu1up, DLU1UP) (const octave_idx_type&, const octave_idx_type&, | |
| 53 double *, const octave_idx_type&, | |
| 54 double *, const octave_idx_type&, | |
| 55 double *, double *); | |
| 56 | |
| 57 F77_RET_T | |
| 58 F77_FUNC (dlup1up, DLUP1UP) (const octave_idx_type&, const octave_idx_type&, | |
| 59 double *, const octave_idx_type&, | |
| 60 double *, const octave_idx_type&, | |
| 61 octave_idx_type *, const double *, | |
| 62 const double *, double *); | |
| 63 #endif | |
| 457 | 64 } |
| 65 | |
| 66 LU::LU (const Matrix& a) | |
| 67 { | |
| 5275 | 68 octave_idx_type a_nr = a.rows (); |
| 69 octave_idx_type a_nc = a.cols (); | |
| 70 octave_idx_type mn = (a_nr < a_nc ? a_nr : a_nc); | |
| 1919 | 71 |
|
10350
12884915a8e4
merge MArray classes & improve Array interface
Jaroslav Hajek <highegg@gmail.com>
parents:
10314
diff
changeset
|
72 ipvt.resize (mn, 1); |
| 5275 | 73 octave_idx_type *pipvt = ipvt.fortran_vec (); |
| 1919 | 74 |
| 1992 | 75 a_fact = a; |
| 76 double *tmp_data = a_fact.fortran_vec (); | |
| 1919 | 77 |
| 5275 | 78 octave_idx_type info = 0; |
| 457 | 79 |
| 4329 | 80 F77_XFCN (dgetrf, DGETRF, (a_nr, a_nc, tmp_data, a_nr, pipvt, info)); |
| 457 | 81 |
|
9694
50db3c5175b5
allow unpacked form of LU
Jaroslav Hajek <highegg@gmail.com>
parents:
8920
diff
changeset
|
82 for (octave_idx_type i = 0; i < mn; i++) |
|
50db3c5175b5
allow unpacked form of LU
Jaroslav Hajek <highegg@gmail.com>
parents:
8920
diff
changeset
|
83 pipvt[i] -= 1; |
| 457 | 84 } |
| 85 | |
| 9708 | 86 #ifdef HAVE_QRUPDATE_LUU |
| 87 | |
| 88 void LU::update (const ColumnVector& u, const ColumnVector& v) | |
| 89 { | |
| 90 if (packed ()) | |
| 91 unpack (); | |
| 92 | |
| 93 Matrix& l = l_fact; | |
| 94 Matrix& r = a_fact; | |
| 95 | |
| 96 octave_idx_type m = l.rows (); | |
| 97 octave_idx_type n = r.columns (); | |
| 98 octave_idx_type k = l.columns (); | |
| 99 | |
| 100 if (u.length () == m && v.length () == n) | |
| 101 { | |
| 102 ColumnVector utmp = u, vtmp = v; | |
| 103 F77_XFCN (dlu1up, DLU1UP, (m, n, l.fortran_vec (), m, r.fortran_vec (), k, | |
| 104 utmp.fortran_vec (), vtmp.fortran_vec ())); | |
| 105 } | |
| 106 else | |
| 107 (*current_liboctave_error_handler) ("luupdate: dimensions mismatch"); | |
| 108 } | |
| 109 | |
| 110 void LU::update (const Matrix& u, const Matrix& v) | |
| 111 { | |
| 112 if (packed ()) | |
| 113 unpack (); | |
| 114 | |
| 115 Matrix& l = l_fact; | |
| 116 Matrix& r = a_fact; | |
| 117 | |
| 118 octave_idx_type m = l.rows (); | |
| 119 octave_idx_type n = r.columns (); | |
| 120 octave_idx_type k = l.columns (); | |
| 121 | |
| 122 if (u.rows () == m && v.rows () == n && u.cols () == v.cols ()) | |
| 123 { | |
| 124 for (volatile octave_idx_type i = 0; i < u.cols (); i++) | |
| 125 { | |
| 126 ColumnVector utmp = u.column (i), vtmp = v.column (i); | |
| 127 F77_XFCN (dlu1up, DLU1UP, (m, n, l.fortran_vec (), m, r.fortran_vec (), k, | |
| 128 utmp.fortran_vec (), vtmp.fortran_vec ())); | |
| 129 } | |
| 130 } | |
| 131 else | |
| 132 (*current_liboctave_error_handler) ("luupdate: dimensions mismatch"); | |
| 133 } | |
| 134 | |
| 135 void LU::update_piv (const ColumnVector& u, const ColumnVector& v) | |
| 136 { | |
| 137 if (packed ()) | |
| 138 unpack (); | |
| 139 | |
| 140 Matrix& l = l_fact; | |
| 141 Matrix& r = a_fact; | |
| 142 | |
| 143 octave_idx_type m = l.rows (); | |
| 144 octave_idx_type n = r.columns (); | |
| 145 octave_idx_type k = l.columns (); | |
| 146 | |
| 147 if (u.length () == m && v.length () == n) | |
| 148 { | |
| 149 ColumnVector utmp = u, vtmp = v; | |
| 150 OCTAVE_LOCAL_BUFFER (double, w, m); | |
| 151 for (octave_idx_type i = 0; i < m; i++) ipvt(i) += 1; // increment | |
| 152 F77_XFCN (dlup1up, DLUP1UP, (m, n, l.fortran_vec (), m, r.fortran_vec (), k, | |
| 153 ipvt.fortran_vec (), utmp.data (), vtmp.data (), w)); | |
| 154 for (octave_idx_type i = 0; i < m; i++) ipvt(i) -= 1; // decrement | |
| 155 } | |
| 156 else | |
| 157 (*current_liboctave_error_handler) ("luupdate: dimensions mismatch"); | |
| 158 } | |
| 159 | |
| 160 void LU::update_piv (const Matrix& u, const Matrix& v) | |
| 161 { | |
| 162 if (packed ()) | |
| 163 unpack (); | |
| 164 | |
| 165 Matrix& l = l_fact; | |
| 166 Matrix& r = a_fact; | |
| 167 | |
| 168 octave_idx_type m = l.rows (); | |
| 169 octave_idx_type n = r.columns (); | |
| 170 octave_idx_type k = l.columns (); | |
| 171 | |
| 172 if (u.rows () == m && v.rows () == n && u.cols () == v.cols ()) | |
| 173 { | |
| 174 OCTAVE_LOCAL_BUFFER (double, w, m); | |
| 175 for (octave_idx_type i = 0; i < m; i++) ipvt(i) += 1; // increment | |
| 176 for (volatile octave_idx_type i = 0; i < u.cols (); i++) | |
| 177 { | |
| 178 ColumnVector utmp = u.column (i), vtmp = v.column (i); | |
| 179 F77_XFCN (dlup1up, DLUP1UP, (m, n, l.fortran_vec (), m, r.fortran_vec (), k, | |
| 180 ipvt.fortran_vec (), utmp.data (), vtmp.data (), w)); | |
| 181 } | |
| 182 for (octave_idx_type i = 0; i < m; i++) ipvt(i) -= 1; // decrement | |
| 183 } | |
| 184 else | |
| 185 (*current_liboctave_error_handler) ("luupdate: dimensions mismatch"); | |
| 186 } | |
| 187 | |
| 188 #else | |
| 189 | |
| 190 void LU::update (const ColumnVector&, const ColumnVector&) | |
| 191 { | |
| 192 (*current_liboctave_error_handler) ("luupdate: not available in this version"); | |
| 193 } | |
| 194 | |
| 195 void LU::update (const Matrix&, const Matrix&) | |
| 196 { | |
| 197 (*current_liboctave_error_handler) ("luupdate: not available in this version"); | |
| 198 } | |
| 199 | |
| 200 void LU::update_piv (const ColumnVector&, const ColumnVector&) | |
| 201 { | |
| 202 (*current_liboctave_error_handler) ("luupdate: not available in this version"); | |
| 203 } | |
| 204 | |
| 205 void LU::update_piv (const Matrix&, const Matrix&) | |
| 206 { | |
| 207 (*current_liboctave_error_handler) ("luupdate: not available in this version"); | |
| 208 } | |
| 209 | |
| 210 #endif |