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