Mercurial > hg > octave-lyh
annotate liboctave/CmplxLU.cc @ 10396:a0b51ac0f88a
optimize accumdim with summation
| author | Jaroslav Hajek <highegg@gmail.com> |
|---|---|
| date | Fri, 05 Mar 2010 12:31:30 +0100 |
| 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 "CmplxLU.h" | |
| 1847 | 30 #include "f77-fcn.h" |
| 457 | 31 #include "lo-error.h" |
| 9708 | 32 #include "oct-locbuf.h" |
| 33 #include "CColVector.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 <ComplexMatrix>; |
| 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 (zgetrf, ZGETRF) (const octave_idx_type&, const octave_idx_type&, Complex*, |
|
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 (zlu1up, ZLU1UP) (const octave_idx_type&, const octave_idx_type&, | |
| 53 Complex *, const octave_idx_type&, | |
| 54 Complex *, const octave_idx_type&, | |
| 55 Complex *, Complex *); | |
| 56 | |
| 57 F77_RET_T | |
| 58 F77_FUNC (zlup1up, ZLUP1UP) (const octave_idx_type&, const octave_idx_type&, | |
| 59 Complex *, const octave_idx_type&, | |
| 60 Complex *, const octave_idx_type&, | |
| 61 octave_idx_type *, const Complex *, | |
| 62 const Complex *, Complex *); | |
| 63 #endif | |
| 457 | 64 } |
| 65 | |
| 66 ComplexLU::ComplexLU (const ComplexMatrix& 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 Complex *tmp_data = a_fact.fortran_vec (); | |
| 1919 | 77 |
| 5275 | 78 octave_idx_type info = 0; |
| 457 | 79 |
| 4329 | 80 F77_XFCN (zgetrf, ZGETRF, (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 ComplexLU::update (const ComplexColumnVector& u, const ComplexColumnVector& v) | |
| 89 { | |
| 90 if (packed ()) | |
| 91 unpack (); | |
| 92 | |
| 93 ComplexMatrix& l = l_fact; | |
| 94 ComplexMatrix& 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 ComplexColumnVector utmp = u, vtmp = v; | |
| 103 F77_XFCN (zlu1up, ZLU1UP, (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 ComplexLU::update (const ComplexMatrix& u, const ComplexMatrix& v) | |
| 111 { | |
| 112 if (packed ()) | |
| 113 unpack (); | |
| 114 | |
| 115 ComplexMatrix& l = l_fact; | |
| 116 ComplexMatrix& 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 ComplexColumnVector utmp = u.column (i), vtmp = v.column (i); | |
| 127 F77_XFCN (zlu1up, ZLU1UP, (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 ComplexLU::update_piv (const ComplexColumnVector& u, const ComplexColumnVector& v) | |
| 136 { | |
| 137 if (packed ()) | |
| 138 unpack (); | |
| 139 | |
| 140 ComplexMatrix& l = l_fact; | |
| 141 ComplexMatrix& 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 ComplexColumnVector utmp = u, vtmp = v; | |
| 150 OCTAVE_LOCAL_BUFFER (Complex, w, m); | |
| 151 for (octave_idx_type i = 0; i < m; i++) ipvt(i) += 1; // increment | |
| 152 F77_XFCN (zlup1up, ZLUP1UP, (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 ComplexLU::update_piv (const ComplexMatrix& u, const ComplexMatrix& v) | |
| 161 { | |
| 162 if (packed ()) | |
| 163 unpack (); | |
| 164 | |
| 165 ComplexMatrix& l = l_fact; | |
| 166 ComplexMatrix& 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 (Complex, 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 ComplexColumnVector utmp = u.column (i), vtmp = v.column (i); | |
| 179 F77_XFCN (zlup1up, ZLUP1UP, (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 ComplexLU::update (const ComplexColumnVector&, const ComplexColumnVector&) | |
| 191 { | |
| 192 (*current_liboctave_error_handler) ("luupdate: not available in this version"); | |
| 193 } | |
| 194 | |
| 195 void ComplexLU::update (const ComplexMatrix&, const ComplexMatrix&) | |
| 196 { | |
| 197 (*current_liboctave_error_handler) ("luupdate: not available in this version"); | |
| 198 } | |
| 199 | |
| 200 void ComplexLU::update_piv (const ComplexColumnVector&, const ComplexColumnVector&) | |
| 201 { | |
| 202 (*current_liboctave_error_handler) ("luupdate: not available in this version"); | |
| 203 } | |
| 204 | |
| 205 void ComplexLU::update_piv (const ComplexMatrix&, const ComplexMatrix&) | |
| 206 { | |
| 207 (*current_liboctave_error_handler) ("luupdate: not available in this version"); | |
| 208 } | |
| 209 | |
| 210 #endif |