Mercurial > octave-nkf
annotate liboctave/dbleLU.cc @ 12312:b10ea6efdc58 release-3-4-x ss-3-3-91
version is now 3.3.91
| author | John W. Eaton <jwe@octave.org> |
|---|---|
| date | Mon, 31 Jan 2011 08:36:58 -0500 |
| parents | a83bad07f7e3 |
| children | 72c96de7a403 |
| rev | line source |
|---|---|
| 457 | 1 /* |
| 2 | |
| 11523 | 3 Copyright (C) 1994-2011 John W. Eaton |
| 9708 | 4 Copyright (C) 2009 VZLU Prague, a.s. |
| 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 | |
| 28 #include "dbleLU.h" | |
| 1847 | 29 #include "f77-fcn.h" |
| 457 | 30 #include "lo-error.h" |
| 9708 | 31 #include "oct-locbuf.h" |
| 32 #include "dColVector.h" | |
| 457 | 33 |
| 1992 | 34 // Instantiate the base LU class for the types we need. |
| 35 | |
| 36 #include <base-lu.h> | |
| 37 #include <base-lu.cc> | |
| 38 | |
|
8367
445d27d79f4e
support permutation matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents:
7482
diff
changeset
|
39 template class base_lu <Matrix>; |
| 1992 | 40 |
| 41 // Define the constructor for this particular derivation. | |
| 42 | |
| 457 | 43 extern "C" |
| 44 { | |
| 4552 | 45 F77_RET_T |
| 11518 | 46 F77_FUNC (dgetrf, DGETRF) (const octave_idx_type&, const octave_idx_type&, |
| 47 double*, const octave_idx_type&, | |
| 48 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 |
|
11574
a83bad07f7e3
attempt better backward compatibility for Array resize functions
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
72 ipvt.resize (dim_vector (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 |