Mercurial > octave

comparison liboctave/numeric/lu.cc @ 31607:aac27ad79be6 stable

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
maint: Re-indent code after switch to using namespace macros. * build-env.h, build-env.in.cc, Cell.h, __betainc__.cc, __eigs__.cc, __ftp__.cc, __ichol__.cc, __ilu__.cc, __isprimelarge__.cc, __magick_read__.cc, __pchip_deriv__.cc, amd.cc, base-text-renderer.cc, base-text-renderer.h, besselj.cc, bitfcns.cc, bsxfun.cc, c-file-ptr-stream.h, call-stack.cc, call-stack.h, ccolamd.cc, cellfun.cc, chol.cc, colamd.cc, dasrt.cc, data.cc, debug.cc, defaults.cc, defaults.h, det.cc, display.cc, display.h, dlmread.cc, dynamic-ld.cc, dynamic-ld.h, ellipj.cc, environment.cc, environment.h, error.cc, error.h, errwarn.h, event-manager.cc, event-manager.h, event-queue.cc, event-queue.h, fcn-info.cc, fcn-info.h, fft.cc, fft2.cc, file-io.cc, filter.cc, find.cc, ft-text-renderer.cc, ft-text-renderer.h, gcd.cc, gl-render.cc, gl-render.h, gl2ps-print.cc, gl2ps-print.h, graphics-toolkit.cc, graphics-toolkit.h, graphics.cc, gsvd.cc, gtk-manager.cc, gtk-manager.h, help.cc, help.h, hook-fcn.cc, hook-fcn.h, input.cc, input.h, interpreter-private.cc, interpreter-private.h, interpreter.cc, interpreter.h, inv.cc, jsondecode.cc, jsonencode.cc, latex-text-renderer.cc, latex-text-renderer.h, load-path.cc, load-path.h, load-save.cc, load-save.h, lookup.cc, ls-hdf5.cc, ls-mat4.cc, ls-mat5.cc, lsode.cc, lu.cc, mappers.cc, matrix_type.cc, max.cc, mex.cc, mexproto.h, mxarray.h, mxtypes.in.h, oct-errno.in.cc, oct-hdf5-types.cc, oct-hist.cc, oct-hist.h, oct-map.cc, oct-map.h, oct-opengl.h, oct-prcstrm.h, oct-process.cc, oct-process.h, oct-stdstrm.h, oct-stream.cc, oct-stream.h, oct-strstrm.h, octave-default-image.h, ordqz.cc, ordschur.cc, pager.cc, pager.h, pinv.cc, pow2.cc, pr-output.cc, psi.cc, qr.cc, quadcc.cc, rand.cc, regexp.cc, settings.cc, settings.h, sighandlers.cc, sighandlers.h, sparse-xpow.cc, sqrtm.cc, stack-frame.cc, stack-frame.h, stream-euler.cc, strfns.cc, svd.cc, syminfo.cc, syminfo.h, symrcm.cc, symrec.cc, symrec.h, symscope.cc, symscope.h, symtab.cc, symtab.h, sysdep.cc, sysdep.h, text-engine.cc, text-engine.h, text-renderer.cc, text-renderer.h, time.cc, toplev.cc, typecast.cc, url-handle-manager.cc, url-handle-manager.h, urlwrite.cc, utils.cc, utils.h, variables.cc, variables.h, xdiv.cc, __delaunayn__.cc, __init_fltk__.cc, __init_gnuplot__.cc, __ode15__.cc, __voronoi__.cc, audioread.cc, convhulln.cc, gzip.cc, cdef-class.cc, cdef-class.h, cdef-fwd.h, cdef-manager.cc, cdef-manager.h, cdef-method.cc, cdef-method.h, cdef-object.cc, cdef-object.h, cdef-package.cc, cdef-package.h, cdef-property.cc, cdef-property.h, cdef-utils.cc, cdef-utils.h, ov-base-diag.cc, ov-base-int.cc, ov-base-mat.cc, ov-base-mat.h, ov-base-scalar.cc, ov-base.cc, ov-base.h, ov-bool-mat.cc, ov-bool-mat.h, ov-bool-sparse.cc, ov-bool.cc, ov-builtin.h, ov-cell.cc, ov-ch-mat.cc, ov-class.cc, ov-class.h, ov-classdef.cc, ov-classdef.h, ov-complex.cc, ov-cx-diag.cc, ov-cx-mat.cc, ov-cx-sparse.cc, ov-dld-fcn.cc, ov-dld-fcn.h, ov-fcn-handle.cc, ov-fcn-handle.h, ov-fcn.h, ov-float.cc, ov-flt-complex.cc, ov-flt-cx-diag.cc, ov-flt-cx-mat.cc, ov-flt-re-diag.cc, ov-flt-re-mat.cc, ov-flt-re-mat.h, ov-intx.h, ov-java.cc, ov-lazy-idx.cc, ov-legacy-range.cc, ov-magic-int.cc, ov-mex-fcn.cc, ov-mex-fcn.h, ov-null-mat.cc, ov-perm.cc, ov-range.cc, ov-re-diag.cc, ov-re-mat.cc, ov-re-mat.h, ov-re-sparse.cc, ov-scalar.cc, ov-str-mat.cc, ov-struct.cc, ov-typeinfo.cc, ov-typeinfo.h, ov-usr-fcn.cc, ov-usr-fcn.h, ov.cc, ov.h, ovl.h, octave.cc, octave.h, op-b-sbm.cc, op-bm-sbm.cc, op-cs-scm.cc, op-fm-fcm.cc, op-fs-fcm.cc, op-s-scm.cc, op-scm-cs.cc, op-scm-s.cc, op-sm-cs.cc, ops.h, anon-fcn-validator.cc, anon-fcn-validator.h, bp-table.cc, bp-table.h, comment-list.cc, comment-list.h, filepos.h, lex.h, oct-lvalue.cc, oct-lvalue.h, parse.h, profiler.cc, profiler.h, pt-anon-scopes.cc, pt-anon-scopes.h, pt-arg-list.cc, pt-arg-list.h, pt-args-block.cc, pt-args-block.h, pt-array-list.cc, pt-array-list.h, pt-assign.cc, pt-assign.h, pt-binop.cc, pt-binop.h, pt-bp.cc, pt-bp.h, pt-cbinop.cc, pt-cbinop.h, pt-cell.cc, pt-cell.h, pt-check.cc, pt-check.h, pt-classdef.cc, pt-classdef.h, pt-cmd.h, pt-colon.cc, pt-colon.h, pt-const.cc, pt-const.h, pt-decl.cc, pt-decl.h, pt-eval.cc, pt-eval.h, pt-except.cc, pt-except.h, pt-exp.cc, pt-exp.h, pt-fcn-handle.cc, pt-fcn-handle.h, pt-id.cc, pt-id.h, pt-idx.cc, pt-idx.h, pt-jump.h, pt-loop.cc, pt-loop.h, pt-mat.cc, pt-mat.h, pt-misc.cc, pt-misc.h, pt-pr-code.cc, pt-pr-code.h, pt-select.cc, pt-select.h, pt-spmd.cc, pt-spmd.h, pt-stmt.cc, pt-stmt.h, pt-tm-const.cc, pt-tm-const.h, pt-unop.cc, pt-unop.h, pt-walk.cc, pt-walk.h, pt.cc, pt.h, token.cc, token.h, Range.cc, Range.h, idx-vector.cc, idx-vector.h, range-fwd.h, CollocWt.cc, CollocWt.h, aepbalance.cc, aepbalance.h, chol.cc, chol.h, gepbalance.cc, gepbalance.h, gsvd.cc, gsvd.h, hess.cc, hess.h, lo-mappers.cc, lo-mappers.h, lo-specfun.cc, lo-specfun.h, lu.cc, lu.h, oct-convn.cc, oct-convn.h, oct-fftw.cc, oct-fftw.h, oct-norm.cc, oct-norm.h, oct-rand.cc, oct-rand.h, oct-spparms.cc, oct-spparms.h, qr.cc, qr.h, qrp.cc, qrp.h, randgamma.cc, randgamma.h, randmtzig.cc, randmtzig.h, randpoisson.cc, randpoisson.h, schur.cc, schur.h, sparse-chol.cc, sparse-chol.h, sparse-lu.cc, sparse-lu.h, sparse-qr.cc, sparse-qr.h, svd.cc, svd.h, child-list.cc, child-list.h, dir-ops.cc, dir-ops.h, file-ops.cc, file-ops.h, file-stat.cc, file-stat.h, lo-sysdep.cc, lo-sysdep.h, lo-sysinfo.cc, lo-sysinfo.h, mach-info.cc, mach-info.h, oct-env.cc, oct-env.h, oct-group.cc, oct-group.h, oct-password.cc, oct-password.h, oct-syscalls.cc, oct-syscalls.h, oct-time.cc, oct-time.h, oct-uname.cc, oct-uname.h, action-container.cc, action-container.h, base-list.h, cmd-edit.cc, cmd-edit.h, cmd-hist.cc, cmd-hist.h, f77-fcn.h, file-info.cc, file-info.h, lo-array-errwarn.cc, lo-array-errwarn.h, lo-hash.cc, lo-hash.h, lo-ieee.h, lo-regexp.cc, lo-regexp.h, lo-utils.cc, lo-utils.h, oct-base64.cc, oct-base64.h, oct-glob.cc, oct-glob.h, oct-inttypes.h, oct-mutex.cc, oct-mutex.h, oct-refcount.h, oct-shlib.cc, oct-shlib.h, oct-sparse.cc, oct-sparse.h, oct-string.h, octave-preserve-stream-state.h, pathsearch.cc, pathsearch.h, quit.cc, quit.h, unwind-prot.cc, unwind-prot.h, url-transfer.cc, url-transfer.h: Re-indent code after switch to using namespace macros.
author Rik <rik@octave.org>
date Thu, 01 Dec 2022 18:02:15 -0800
parents e88a07dec498
children 597f3ee61a48
comparison
equal deleted inserted replaced
31605:e88a07dec498 31607:aac27ad79be6
46 46
47 OCTAVE_BEGIN_NAMESPACE(octave) 47 OCTAVE_BEGIN_NAMESPACE(octave)
48 48
49 OCTAVE_BEGIN_NAMESPACE(math) 49 OCTAVE_BEGIN_NAMESPACE(math)
50 50
51 // FIXME: PermMatrix::col_perm_vec returns Array<octave_idx_type> 51 // FIXME: PermMatrix::col_perm_vec returns Array<octave_idx_type>
52 // but m_ipvt is an Array<octave_f77_int_type>. This could cause 52 // but m_ipvt is an Array<octave_f77_int_type>. This could cause
53 // trouble for large arrays if octave_f77_int_type is 32-bits but 53 // trouble for large arrays if octave_f77_int_type is 32-bits but
54 // octave_idx_type is 64. Since this constructor is called from 54 // octave_idx_type is 64. Since this constructor is called from
55 // Fluupdate, it could be given values that are out of range. We 55 // Fluupdate, it could be given values that are out of range. We
56 // should ensure that the values are within range here. 56 // should ensure that the values are within range here.
57 57
58 template <typename T> 58 template <typename T>
59 lu<T>::lu (const T& l, const T& u, const PermMatrix& p) 59 lu<T>::lu (const T& l, const T& u, const PermMatrix& p)
60 : m_a_fact (u), m_L (l), m_ipvt (p.transpose ().col_perm_vec ()) 60 : m_a_fact (u), m_L (l), m_ipvt (p.transpose ().col_perm_vec ())
61 { 61 {
62 if (l.columns () != u.rows ()) 62 if (l.columns () != u.rows ())
63 (*current_liboctave_error_handler) ("lu: dimension mismatch"); 63 (*current_liboctave_error_handler) ("lu: dimension mismatch");
64 } 64 }
65 65
66 template <typename T> 66 template <typename T>
67 bool 67 bool
68 lu<T>::packed (void) const 68 lu<T>::packed (void) const
69 { 69 {
70 return m_L.dims () == dim_vector (); 70 return m_L.dims () == dim_vector ();
71 } 71 }
72 72
73 template <typename T> 73 template <typename T>
74 void 74 void
75 lu<T>::unpack (void) 75 lu<T>::unpack (void)
76 { 76 {
77 if (packed ()) 77 if (packed ())
78 {
79 m_L = L ();
80 m_a_fact = U (); // FIXME: sub-optimal
81
82 // FIXME: getp returns Array<octave_idx_type> but m_ipvt is
83 // Array<octave_f77_int_type>. However, getp produces its
84 // result from a valid m_ipvt array so validation should not be
85 // necessary. OTOH, it might be better to have a version of
86 // getp that doesn't cause us to convert from
87 // Array<octave_f77_int_type> to Array<octave_idx_type> and
88 // back again.
89
90 m_ipvt = getp ();
91 }
92 }
93
94 template <typename T>
95 T
96 lu<T>::L (void) const
97 {
98 if (packed ())
99 {
100 octave_idx_type a_nr = m_a_fact.rows ();
101 octave_idx_type a_nc = m_a_fact.columns ();
102 octave_idx_type mn = (a_nr < a_nc ? a_nr : a_nc);
103
104 T l (a_nr, mn, ELT_T (0.0));
105
106 for (octave_idx_type i = 0; i < a_nr; i++)
78 { 107 {
79 m_L = L (); 108 if (i < a_nc)
80 m_a_fact = U (); // FIXME: sub-optimal 109 l.xelem (i, i) = 1.0;
81 110
82 // FIXME: getp returns Array<octave_idx_type> but m_ipvt is 111 for (octave_idx_type j = 0; j < (i < a_nc ? i : a_nc); j++)
83 // Array<octave_f77_int_type>. However, getp produces its 112 l.xelem (i, j) = m_a_fact.xelem (i, j);
84 // result from a valid m_ipvt array so validation should not be
85 // necessary. OTOH, it might be better to have a version of
86 // getp that doesn't cause us to convert from
87 // Array<octave_f77_int_type> to Array<octave_idx_type> and
88 // back again.
89
90 m_ipvt = getp ();
91 } 113 }
92 } 114
93 115 return l;
94 template <typename T> 116 }
95 T 117 else
96 lu<T>::L (void) const 118 return m_L;
97 { 119 }
98 if (packed ()) 120
121 template <typename T>
122 T
123 lu<T>::U (void) const
124 {
125 if (packed ())
126 {
127 octave_idx_type a_nr = m_a_fact.rows ();
128 octave_idx_type a_nc = m_a_fact.columns ();
129 octave_idx_type mn = (a_nr < a_nc ? a_nr : a_nc);
130
131 T u (mn, a_nc, ELT_T (0.0));
132
133 for (octave_idx_type i = 0; i < mn; i++)
99 { 134 {
100 octave_idx_type a_nr = m_a_fact.rows (); 135 for (octave_idx_type j = i; j < a_nc; j++)
101 octave_idx_type a_nc = m_a_fact.columns (); 136 u.xelem (i, j) = m_a_fact.xelem (i, j);
102 octave_idx_type mn = (a_nr < a_nc ? a_nr : a_nc); 137 }
103 138
104 T l (a_nr, mn, ELT_T (0.0)); 139 return u;
105 140 }
106 for (octave_idx_type i = 0; i < a_nr; i++) 141 else
142 return m_a_fact;
143 }
144
145 template <typename T>
146 T
147 lu<T>::Y (void) const
148 {
149 if (! packed ())
150 (*current_liboctave_error_handler)
151 ("lu: Y () not implemented for unpacked form");
152
153 return m_a_fact;
154 }
155
156 template <typename T>
157 Array<octave_idx_type>
158 lu<T>::getp (void) const
159 {
160 if (packed ())
161 {
162 octave_idx_type a_nr = m_a_fact.rows ();
163
164 Array<octave_idx_type> pvt (dim_vector (a_nr, 1));
165
166 for (octave_idx_type i = 0; i < a_nr; i++)
167 pvt.xelem (i) = i;
168
169 for (octave_idx_type i = 0; i < m_ipvt.numel (); i++)
170 {
171 octave_idx_type k = m_ipvt.xelem (i);
172
173 if (k != i)
107 { 174 {
108 if (i < a_nc) 175 octave_idx_type tmp = pvt.xelem (k);
109 l.xelem (i, i) = 1.0; 176 pvt.xelem (k) = pvt.xelem (i);
110 177 pvt.xelem (i) = tmp;
111 for (octave_idx_type j = 0; j < (i < a_nc ? i : a_nc); j++)
112 l.xelem (i, j) = m_a_fact.xelem (i, j);
113 }
114
115 return l;
116 }
117 else
118 return m_L;
119 }
120
121 template <typename T>
122 T
123 lu<T>::U (void) const
124 {
125 if (packed ())
126 {
127 octave_idx_type a_nr = m_a_fact.rows ();
128 octave_idx_type a_nc = m_a_fact.columns ();
129 octave_idx_type mn = (a_nr < a_nc ? a_nr : a_nc);
130
131 T u (mn, a_nc, ELT_T (0.0));
132
133 for (octave_idx_type i = 0; i < mn; i++)
134 {
135 for (octave_idx_type j = i; j < a_nc; j++)
136 u.xelem (i, j) = m_a_fact.xelem (i, j);
137 }
138
139 return u;
140 }
141 else
142 return m_a_fact;
143 }
144
145 template <typename T>
146 T
147 lu<T>::Y (void) const
148 {
149 if (! packed ())
150 (*current_liboctave_error_handler)
151 ("lu: Y () not implemented for unpacked form");
152
153 return m_a_fact;
154 }
155
156 template <typename T>
157 Array<octave_idx_type>
158 lu<T>::getp (void) const
159 {
160 if (packed ())
161 {
162 octave_idx_type a_nr = m_a_fact.rows ();
163
164 Array<octave_idx_type> pvt (dim_vector (a_nr, 1));
165
166 for (octave_idx_type i = 0; i < a_nr; i++)
167 pvt.xelem (i) = i;
168
169 for (octave_idx_type i = 0; i < m_ipvt.numel (); i++)
170 {
171 octave_idx_type k = m_ipvt.xelem (i);
172
173 if (k != i)
174 {
175 octave_idx_type tmp = pvt.xelem (k);
176 pvt.xelem (k) = pvt.xelem (i);
177 pvt.xelem (i) = tmp;
178 }
179 }
180
181 return pvt;
182 }
183 else
184 return m_ipvt;
185 }
186
187 template <typename T>
188 PermMatrix
189 lu<T>::P (void) const
190 {
191 return PermMatrix (getp (), false);
192 }
193
194 template <typename T>
195 ColumnVector
196 lu<T>::P_vec (void) const
197 {
198 octave_idx_type a_nr = m_a_fact.rows ();
199
200 ColumnVector p (a_nr);
201
202 Array<octave_idx_type> pvt = getp ();
203
204 for (octave_idx_type i = 0; i < a_nr; i++)
205 p.xelem (i) = static_cast<double> (pvt.xelem (i) + 1);
206
207 return p;
208 }
209
210 template <typename T>
211 bool
212 lu<T>::regular (void) const
213 {
214 bool retval = true;
215
216 octave_idx_type k = std::min (m_a_fact.rows (), m_a_fact.columns ());
217
218 for (octave_idx_type i = 0; i < k; i++)
219 {
220 if (m_a_fact(i, i) == ELT_T ())
221 {
222 retval = false;
223 break;
224 } 178 }
225 } 179 }
226 180
227 return retval; 181 return pvt;
228 } 182 }
183 else
184 return m_ipvt;
185 }
186
187 template <typename T>
188 PermMatrix
189 lu<T>::P (void) const
190 {
191 return PermMatrix (getp (), false);
192 }
193
194 template <typename T>
195 ColumnVector
196 lu<T>::P_vec (void) const
197 {
198 octave_idx_type a_nr = m_a_fact.rows ();
199
200 ColumnVector p (a_nr);
201
202 Array<octave_idx_type> pvt = getp ();
203
204 for (octave_idx_type i = 0; i < a_nr; i++)
205 p.xelem (i) = static_cast<double> (pvt.xelem (i) + 1);
206
207 return p;
208 }
209
210 template <typename T>
211 bool
212 lu<T>::regular (void) const
213 {
214 bool retval = true;
215
216 octave_idx_type k = std::min (m_a_fact.rows (), m_a_fact.columns ());
217
218 for (octave_idx_type i = 0; i < k; i++)
219 {
220 if (m_a_fact(i, i) == ELT_T ())
221 {
222 retval = false;
223 break;
224 }
225 }
226
227 return retval;
228 }
229 229
230 #if ! defined (HAVE_QRUPDATE_LUU) 230 #if ! defined (HAVE_QRUPDATE_LUU)
231 231
232 template <typename T> 232 template <typename T>
233 void 233 void
234 lu<T>::update (const VT&, const VT&) 234 lu<T>::update (const VT&, const VT&)
235 { 235 {
236 (*current_liboctave_error_handler) 236 (*current_liboctave_error_handler)
237 ("luupdate: support for qrupdate with LU updates " 237 ("luupdate: support for qrupdate with LU updates "
238 "was unavailable or disabled when liboctave was built"); 238 "was unavailable or disabled when liboctave was built");
239 } 239 }
240 240
241 template <typename T> 241 template <typename T>
242 void 242 void
243 lu<T>::update (const T&, const T&) 243 lu<T>::update (const T&, const T&)
244 { 244 {
245 (*current_liboctave_error_handler) 245 (*current_liboctave_error_handler)
246 ("luupdate: support for qrupdate with LU updates " 246 ("luupdate: support for qrupdate with LU updates "
247 "was unavailable or disabled when liboctave was built"); 247 "was unavailable or disabled when liboctave was built");
248 } 248 }
249 249
250 template <typename T> 250 template <typename T>
251 void 251 void
252 lu<T>::update_piv (const VT&, const VT&) 252 lu<T>::update_piv (const VT&, const VT&)
253 { 253 {
254 (*current_liboctave_error_handler) 254 (*current_liboctave_error_handler)
255 ("luupdate: support for qrupdate with LU updates " 255 ("luupdate: support for qrupdate with LU updates "
256 "was unavailable or disabled when liboctave was built"); 256 "was unavailable or disabled when liboctave was built");
257 } 257 }
258 258
259 template <typename T> 259 template <typename T>
260 void 260 void
261 lu<T>::update_piv (const T&, const T&) 261 lu<T>::update_piv (const T&, const T&)
262 { 262 {
263 (*current_liboctave_error_handler) 263 (*current_liboctave_error_handler)
264 ("luupdate: support for qrupdate with LU updates " 264 ("luupdate: support for qrupdate with LU updates "
265 "was unavailable or disabled when liboctave was built"); 265 "was unavailable or disabled when liboctave was built");
266 } 266 }
267 267
268 #endif 268 #endif
269 269
270 // Specializations. 270 // Specializations.
271 271
272 template <> 272 template <>
273 OCTAVE_API 273 OCTAVE_API
274 lu<Matrix>::lu (const Matrix& a) 274 lu<Matrix>::lu (const Matrix& a)
275 { 275 {
276 F77_INT a_nr = to_f77_int (a.rows ()); 276 F77_INT a_nr = to_f77_int (a.rows ());
277 F77_INT a_nc = to_f77_int (a.columns ()); 277 F77_INT a_nc = to_f77_int (a.columns ());
278 F77_INT mn = (a_nr < a_nc ? a_nr : a_nc); 278 F77_INT mn = (a_nr < a_nc ? a_nr : a_nc);
279 279
280 m_ipvt.resize (dim_vector (mn, 1)); 280 m_ipvt.resize (dim_vector (mn, 1));
281 F77_INT *pipvt = m_ipvt.fortran_vec (); 281 F77_INT *pipvt = m_ipvt.fortran_vec ();
282 282
283 m_a_fact = a; 283 m_a_fact = a;
284 double *tmp_data = m_a_fact.fortran_vec (); 284 double *tmp_data = m_a_fact.fortran_vec ();
285 285
286 F77_INT info = 0; 286 F77_INT info = 0;
287 287
288 F77_XFCN (dgetrf, DGETRF, (a_nr, a_nc, tmp_data, a_nr, pipvt, info)); 288 F77_XFCN (dgetrf, DGETRF, (a_nr, a_nc, tmp_data, a_nr, pipvt, info));
289 289
290 for (F77_INT i = 0; i < mn; i++) 290 for (F77_INT i = 0; i < mn; i++)
291 pipvt[i] -= 1; 291 pipvt[i] -= 1;
292 } 292 }
293 293
294 #if defined (HAVE_QRUPDATE_LUU) 294 #if defined (HAVE_QRUPDATE_LUU)
295 295
296 template <> 296 template <>
297 OCTAVE_API void 297 OCTAVE_API void
298 lu<Matrix>::update (const ColumnVector& u, const ColumnVector& v) 298 lu<Matrix>::update (const ColumnVector& u, const ColumnVector& v)
299 { 299 {
300 if (packed ()) 300 if (packed ())
301 unpack (); 301 unpack ();
302 302
303 Matrix& l = m_L; 303 Matrix& l = m_L;
304 Matrix& r = m_a_fact; 304 Matrix& r = m_a_fact;
305 305
306 F77_INT m = to_f77_int (l.rows ()); 306 F77_INT m = to_f77_int (l.rows ());
307 F77_INT n = to_f77_int (r.columns ()); 307 F77_INT n = to_f77_int (r.columns ());
308 F77_INT k = to_f77_int (l.columns ()); 308 F77_INT k = to_f77_int (l.columns ());
309 309
310 F77_INT u_nel = to_f77_int (u.numel ()); 310 F77_INT u_nel = to_f77_int (u.numel ());
311 F77_INT v_nel = to_f77_int (v.numel ()); 311 F77_INT v_nel = to_f77_int (v.numel ());
312 312
313 if (u_nel != m || v_nel != n) 313 if (u_nel != m || v_nel != n)
314 (*current_liboctave_error_handler) ("luupdate: dimensions mismatch"); 314 (*current_liboctave_error_handler) ("luupdate: dimensions mismatch");
315 315
316 ColumnVector utmp = u; 316 ColumnVector utmp = u;
317 ColumnVector vtmp = v; 317 ColumnVector vtmp = v;
318 F77_XFCN (dlu1up, DLU1UP, (m, n, l.fortran_vec (), m, r.fortran_vec (), 318 F77_XFCN (dlu1up, DLU1UP, (m, n, l.fortran_vec (), m, r.fortran_vec (),
319 k, utmp.fortran_vec (), vtmp.fortran_vec ())); 319 k, utmp.fortran_vec (), vtmp.fortran_vec ()));
320 } 320 }
321 321
322 template <> 322 template <>
323 OCTAVE_API void 323 OCTAVE_API void
324 lu<Matrix>::update (const Matrix& u, const Matrix& v) 324 lu<Matrix>::update (const Matrix& u, const Matrix& v)
325 { 325 {
326 if (packed ()) 326 if (packed ())
327 unpack (); 327 unpack ();
328 328
329 Matrix& l = m_L; 329 Matrix& l = m_L;
330 Matrix& r = m_a_fact; 330 Matrix& r = m_a_fact;
331 331
332 F77_INT m = to_f77_int (l.rows ()); 332 F77_INT m = to_f77_int (l.rows ());
333 F77_INT n = to_f77_int (r.columns ()); 333 F77_INT n = to_f77_int (r.columns ());
334 F77_INT k = to_f77_int (l.columns ()); 334 F77_INT k = to_f77_int (l.columns ());
335 335
336 F77_INT u_nr = to_f77_int (u.rows ()); 336 F77_INT u_nr = to_f77_int (u.rows ());
337 F77_INT u_nc = to_f77_int (u.columns ()); 337 F77_INT u_nc = to_f77_int (u.columns ());
338 338
339 F77_INT v_nr = to_f77_int (v.rows ()); 339 F77_INT v_nr = to_f77_int (v.rows ());
340 F77_INT v_nc = to_f77_int (v.columns ()); 340 F77_INT v_nc = to_f77_int (v.columns ());
341 341
342 if (u_nr != m || v_nr != n || u_nc != v_nc) 342 if (u_nr != m || v_nr != n || u_nc != v_nc)
343 (*current_liboctave_error_handler) ("luupdate: dimensions mismatch"); 343 (*current_liboctave_error_handler) ("luupdate: dimensions mismatch");
344 344
345 for (volatile F77_INT i = 0; i < u_nc; i++) 345 for (volatile F77_INT i = 0; i < u_nc; i++)
346 { 346 {
347 ColumnVector utmp = u.column (i); 347 ColumnVector utmp = u.column (i);
348 ColumnVector vtmp = v.column (i); 348 ColumnVector vtmp = v.column (i);
349 F77_XFCN (dlu1up, DLU1UP, (m, n, l.fortran_vec (), 349 F77_XFCN (dlu1up, DLU1UP, (m, n, l.fortran_vec (),
350 m, r.fortran_vec (), k, 350 m, r.fortran_vec (), k,
351 utmp.fortran_vec (), vtmp.fortran_vec ())); 351 utmp.fortran_vec (), vtmp.fortran_vec ()));
352 } 352 }
353 } 353 }
354 354
355 template <> 355 template <>
356 OCTAVE_API void 356 OCTAVE_API void
357 lu<Matrix>::update_piv (const ColumnVector& u, const ColumnVector& v) 357 lu<Matrix>::update_piv (const ColumnVector& u, const ColumnVector& v)
358 { 358 {
359 if (packed ()) 359 if (packed ())
360 unpack (); 360 unpack ();
361 361
362 Matrix& l = m_L; 362 Matrix& l = m_L;
363 Matrix& r = m_a_fact; 363 Matrix& r = m_a_fact;
364 364
365 F77_INT m = to_f77_int (l.rows ()); 365 F77_INT m = to_f77_int (l.rows ());
366 F77_INT n = to_f77_int (r.columns ()); 366 F77_INT n = to_f77_int (r.columns ());
367 F77_INT k = to_f77_int (l.columns ()); 367 F77_INT k = to_f77_int (l.columns ());
368 368
369 F77_INT u_nel = to_f77_int (u.numel ()); 369 F77_INT u_nel = to_f77_int (u.numel ());
370 F77_INT v_nel = to_f77_int (v.numel ()); 370 F77_INT v_nel = to_f77_int (v.numel ());
371 371
372 if (u_nel != m || v_nel != n) 372 if (u_nel != m || v_nel != n)
373 (*current_liboctave_error_handler) ("luupdate: dimensions mismatch"); 373 (*current_liboctave_error_handler) ("luupdate: dimensions mismatch");
374 374
375 ColumnVector utmp = u; 375 ColumnVector utmp = u;
376 ColumnVector vtmp = v; 376 ColumnVector vtmp = v;
377 OCTAVE_LOCAL_BUFFER (double, w, m); 377 OCTAVE_LOCAL_BUFFER (double, w, m);
378 for (F77_INT i = 0; i < m; i++) m_ipvt(i) += 1; // increment 378 for (F77_INT i = 0; i < m; i++) m_ipvt(i) += 1; // increment
379 F77_XFCN (dlup1up, DLUP1UP, (m, n, l.fortran_vec (),
380 m, r.fortran_vec (), k,
381 m_ipvt.fortran_vec (),
382 utmp.data (), vtmp.data (), w));
383 for (F77_INT i = 0; i < m; i++) m_ipvt(i) -= 1; // decrement
384 }
385
386 template <>
387 OCTAVE_API void
388 lu<Matrix>::update_piv (const Matrix& u, const Matrix& v)
389 {
390 if (packed ())
391 unpack ();
392
393 Matrix& l = m_L;
394 Matrix& r = m_a_fact;
395
396 F77_INT m = to_f77_int (l.rows ());
397 F77_INT n = to_f77_int (r.columns ());
398 F77_INT k = to_f77_int (l.columns ());
399
400 F77_INT u_nr = to_f77_int (u.rows ());
401 F77_INT u_nc = to_f77_int (u.columns ());
402
403 F77_INT v_nr = to_f77_int (v.rows ());
404 F77_INT v_nc = to_f77_int (v.columns ());
405
406 if (u_nr != m || v_nr != n || u_nc != v_nc)
407 (*current_liboctave_error_handler) ("luupdate: dimensions mismatch");
408
409 OCTAVE_LOCAL_BUFFER (double, w, m);
410 for (F77_INT i = 0; i < m; i++) m_ipvt(i) += 1; // increment
411 for (volatile F77_INT i = 0; i < u_nc; i++)
412 {
413 ColumnVector utmp = u.column (i);
414 ColumnVector vtmp = v.column (i);
379 F77_XFCN (dlup1up, DLUP1UP, (m, n, l.fortran_vec (), 415 F77_XFCN (dlup1up, DLUP1UP, (m, n, l.fortran_vec (),
380 m, r.fortran_vec (), k, 416 m, r.fortran_vec (), k,
381 m_ipvt.fortran_vec (), 417 m_ipvt.fortran_vec (),
382 utmp.data (), vtmp.data (), w)); 418 utmp.data (), vtmp.data (), w));
383 for (F77_INT i = 0; i < m; i++) m_ipvt(i) -= 1; // decrement 419 }
384 } 420 for (F77_INT i = 0; i < m; i++) m_ipvt(i) -= 1; // decrement
385 421 }
386 template <>
387 OCTAVE_API void
388 lu<Matrix>::update_piv (const Matrix& u, const Matrix& v)
389 {
390 if (packed ())
391 unpack ();
392
393 Matrix& l = m_L;
394 Matrix& r = m_a_fact;
395
396 F77_INT m = to_f77_int (l.rows ());
397 F77_INT n = to_f77_int (r.columns ());
398 F77_INT k = to_f77_int (l.columns ());
399
400 F77_INT u_nr = to_f77_int (u.rows ());
401 F77_INT u_nc = to_f77_int (u.columns ());
402
403 F77_INT v_nr = to_f77_int (v.rows ());
404 F77_INT v_nc = to_f77_int (v.columns ());
405
406 if (u_nr != m || v_nr != n || u_nc != v_nc)
407 (*current_liboctave_error_handler) ("luupdate: dimensions mismatch");
408
409 OCTAVE_LOCAL_BUFFER (double, w, m);
410 for (F77_INT i = 0; i < m; i++) m_ipvt(i) += 1; // increment
411 for (volatile F77_INT i = 0; i < u_nc; i++)
412 {
413 ColumnVector utmp = u.column (i);
414 ColumnVector vtmp = v.column (i);
415 F77_XFCN (dlup1up, DLUP1UP, (m, n, l.fortran_vec (),
416 m, r.fortran_vec (), k,
417 m_ipvt.fortran_vec (),
418 utmp.data (), vtmp.data (), w));
419 }
420 for (F77_INT i = 0; i < m; i++) m_ipvt(i) -= 1; // decrement
421 }
422 422
423 #endif 423 #endif
424 424
425 template <> 425 template <>
426 OCTAVE_API 426 OCTAVE_API
427 lu<FloatMatrix>::lu (const FloatMatrix& a) 427 lu<FloatMatrix>::lu (const FloatMatrix& a)
428 { 428 {
429 F77_INT a_nr = to_f77_int (a.rows ()); 429 F77_INT a_nr = to_f77_int (a.rows ());
430 F77_INT a_nc = to_f77_int (a.columns ()); 430 F77_INT a_nc = to_f77_int (a.columns ());
431 F77_INT mn = (a_nr < a_nc ? a_nr : a_nc); 431 F77_INT mn = (a_nr < a_nc ? a_nr : a_nc);
432 432
433 m_ipvt.resize (dim_vector (mn, 1)); 433 m_ipvt.resize (dim_vector (mn, 1));
434 F77_INT *pipvt = m_ipvt.fortran_vec (); 434 F77_INT *pipvt = m_ipvt.fortran_vec ();
435 435
436 m_a_fact = a; 436 m_a_fact = a;
437 float *tmp_data = m_a_fact.fortran_vec (); 437 float *tmp_data = m_a_fact.fortran_vec ();
438 438
439 F77_INT info = 0; 439 F77_INT info = 0;
440 440
441 F77_XFCN (sgetrf, SGETRF, (a_nr, a_nc, tmp_data, a_nr, pipvt, info)); 441 F77_XFCN (sgetrf, SGETRF, (a_nr, a_nc, tmp_data, a_nr, pipvt, info));
442 442
443 for (F77_INT i = 0; i < mn; i++) 443 for (F77_INT i = 0; i < mn; i++)
444 pipvt[i] -= 1; 444 pipvt[i] -= 1;
445 } 445 }
446 446
447 #if defined (HAVE_QRUPDATE_LUU) 447 #if defined (HAVE_QRUPDATE_LUU)
448 448
449 template <> 449 template <>
450 OCTAVE_API void 450 OCTAVE_API void
451 lu<FloatMatrix>::update (const FloatColumnVector& u, 451 lu<FloatMatrix>::update (const FloatColumnVector& u,
452 const FloatColumnVector& v) 452 const FloatColumnVector& v)
453 { 453 {
454 if (packed ()) 454 if (packed ())
455 unpack (); 455 unpack ();
456 456
457 FloatMatrix& l = m_L; 457 FloatMatrix& l = m_L;
458 FloatMatrix& r = m_a_fact; 458 FloatMatrix& r = m_a_fact;
459 459
460 F77_INT m = to_f77_int (l.rows ()); 460 F77_INT m = to_f77_int (l.rows ());
461 F77_INT n = to_f77_int (r.columns ()); 461 F77_INT n = to_f77_int (r.columns ());
462 F77_INT k = to_f77_int (l.columns ()); 462 F77_INT k = to_f77_int (l.columns ());
463 463
464 F77_INT u_nel = to_f77_int (u.numel ()); 464 F77_INT u_nel = to_f77_int (u.numel ());
465 F77_INT v_nel = to_f77_int (v.numel ()); 465 F77_INT v_nel = to_f77_int (v.numel ());
466 466
467 if (u_nel != m || v_nel != n) 467 if (u_nel != m || v_nel != n)
468 (*current_liboctave_error_handler) ("luupdate: dimensions mismatch"); 468 (*current_liboctave_error_handler) ("luupdate: dimensions mismatch");
469 469
470 FloatColumnVector utmp = u; 470 FloatColumnVector utmp = u;
471 FloatColumnVector vtmp = v; 471 FloatColumnVector vtmp = v;
472 F77_XFCN (slu1up, SLU1UP, (m, n, l.fortran_vec (),
473 m, r.fortran_vec (), k,
474 utmp.fortran_vec (), vtmp.fortran_vec ()));
475 }
476
477 template <>
478 OCTAVE_API void
479 lu<FloatMatrix>::update (const FloatMatrix& u, const FloatMatrix& v)
480 {
481 if (packed ())
482 unpack ();
483
484 FloatMatrix& l = m_L;
485 FloatMatrix& r = m_a_fact;
486
487 F77_INT m = to_f77_int (l.rows ());
488 F77_INT n = to_f77_int (r.columns ());
489 F77_INT k = to_f77_int (l.columns ());
490
491 F77_INT u_nr = to_f77_int (u.rows ());
492 F77_INT u_nc = to_f77_int (u.columns ());
493
494 F77_INT v_nr = to_f77_int (v.rows ());
495 F77_INT v_nc = to_f77_int (v.columns ());
496
497 if (u_nr != m || v_nr != n || u_nc != v_nc)
498 (*current_liboctave_error_handler) ("luupdate: dimensions mismatch");
499
500 for (volatile F77_INT i = 0; i < u_nc; i++)
501 {
502 FloatColumnVector utmp = u.column (i);
503 FloatColumnVector vtmp = v.column (i);
472 F77_XFCN (slu1up, SLU1UP, (m, n, l.fortran_vec (), 504 F77_XFCN (slu1up, SLU1UP, (m, n, l.fortran_vec (),
473 m, r.fortran_vec (), k, 505 m, r.fortran_vec (), k,
474 utmp.fortran_vec (), vtmp.fortran_vec ())); 506 utmp.fortran_vec (), vtmp.fortran_vec ()));
475 } 507 }
476 508 }
477 template <> 509
478 OCTAVE_API void 510 template <>
479 lu<FloatMatrix>::update (const FloatMatrix& u, const FloatMatrix& v) 511 OCTAVE_API void
480 { 512 lu<FloatMatrix>::update_piv (const FloatColumnVector& u,
481 if (packed ()) 513 const FloatColumnVector& v)
482 unpack (); 514 {
483 515 if (packed ())
484 FloatMatrix& l = m_L; 516 unpack ();
485 FloatMatrix& r = m_a_fact; 517
486 518 FloatMatrix& l = m_L;
487 F77_INT m = to_f77_int (l.rows ()); 519 FloatMatrix& r = m_a_fact;
488 F77_INT n = to_f77_int (r.columns ()); 520
489 F77_INT k = to_f77_int (l.columns ()); 521 F77_INT m = to_f77_int (l.rows ());
490 522 F77_INT n = to_f77_int (r.columns ());
491 F77_INT u_nr = to_f77_int (u.rows ()); 523 F77_INT k = to_f77_int (l.columns ());
492 F77_INT u_nc = to_f77_int (u.columns ()); 524
493 525 F77_INT u_nel = to_f77_int (u.numel ());
494 F77_INT v_nr = to_f77_int (v.rows ()); 526 F77_INT v_nel = to_f77_int (v.numel ());
495 F77_INT v_nc = to_f77_int (v.columns ()); 527
496 528 if (u_nel != m || v_nel != n)
497 if (u_nr != m || v_nr != n || u_nc != v_nc) 529 (*current_liboctave_error_handler) ("luupdate: dimensions mismatch");
498 (*current_liboctave_error_handler) ("luupdate: dimensions mismatch"); 530
499 531 FloatColumnVector utmp = u;
500 for (volatile F77_INT i = 0; i < u_nc; i++) 532 FloatColumnVector vtmp = v;
501 { 533 OCTAVE_LOCAL_BUFFER (float, w, m);
502 FloatColumnVector utmp = u.column (i); 534 for (F77_INT i = 0; i < m; i++) m_ipvt(i) += 1; // increment
503 FloatColumnVector vtmp = v.column (i); 535 F77_XFCN (slup1up, SLUP1UP, (m, n, l.fortran_vec (),
504 F77_XFCN (slu1up, SLU1UP, (m, n, l.fortran_vec (), 536 m, r.fortran_vec (), k,
505 m, r.fortran_vec (), k, 537 m_ipvt.fortran_vec (),
506 utmp.fortran_vec (), vtmp.fortran_vec ())); 538 utmp.data (), vtmp.data (), w));
507 } 539 for (F77_INT i = 0; i < m; i++) m_ipvt(i) -= 1; // decrement
508 } 540 }
509 541
510 template <> 542 template <>
511 OCTAVE_API void 543 OCTAVE_API void
512 lu<FloatMatrix>::update_piv (const FloatColumnVector& u, 544 lu<FloatMatrix>::update_piv (const FloatMatrix& u, const FloatMatrix& v)
513 const FloatColumnVector& v) 545 {
514 { 546 if (packed ())
515 if (packed ()) 547 unpack ();
516 unpack (); 548
517 549 FloatMatrix& l = m_L;
518 FloatMatrix& l = m_L; 550 FloatMatrix& r = m_a_fact;
519 FloatMatrix& r = m_a_fact; 551
520 552 F77_INT m = to_f77_int (l.rows ());
521 F77_INT m = to_f77_int (l.rows ()); 553 F77_INT n = to_f77_int (r.columns ());
522 F77_INT n = to_f77_int (r.columns ()); 554 F77_INT k = to_f77_int (l.columns ());
523 F77_INT k = to_f77_int (l.columns ()); 555
524 556 F77_INT u_nr = to_f77_int (u.rows ());
525 F77_INT u_nel = to_f77_int (u.numel ()); 557 F77_INT u_nc = to_f77_int (u.columns ());
526 F77_INT v_nel = to_f77_int (v.numel ()); 558
527 559 F77_INT v_nr = to_f77_int (v.rows ());
528 if (u_nel != m || v_nel != n) 560 F77_INT v_nc = to_f77_int (v.columns ());
529 (*current_liboctave_error_handler) ("luupdate: dimensions mismatch"); 561
530 562 if (u_nr != m || v_nr != n || u_nc != v_nc)
531 FloatColumnVector utmp = u; 563 (*current_liboctave_error_handler) ("luupdate: dimensions mismatch");
532 FloatColumnVector vtmp = v; 564
533 OCTAVE_LOCAL_BUFFER (float, w, m); 565 OCTAVE_LOCAL_BUFFER (float, w, m);
534 for (F77_INT i = 0; i < m; i++) m_ipvt(i) += 1; // increment 566 for (F77_INT i = 0; i < m; i++) m_ipvt(i) += 1; // increment
567 for (volatile F77_INT i = 0; i < u_nc; i++)
568 {
569 FloatColumnVector utmp = u.column (i);
570 FloatColumnVector vtmp = v.column (i);
535 F77_XFCN (slup1up, SLUP1UP, (m, n, l.fortran_vec (), 571 F77_XFCN (slup1up, SLUP1UP, (m, n, l.fortran_vec (),
536 m, r.fortran_vec (), k, 572 m, r.fortran_vec (), k,
537 m_ipvt.fortran_vec (), 573 m_ipvt.fortran_vec (),
538 utmp.data (), vtmp.data (), w)); 574 utmp.data (), vtmp.data (), w));
539 for (F77_INT i = 0; i < m; i++) m_ipvt(i) -= 1; // decrement 575 }
540 } 576 for (F77_INT i = 0; i < m; i++) m_ipvt(i) -= 1; // decrement
541 577 }
542 template <>
543 OCTAVE_API void
544 lu<FloatMatrix>::update_piv (const FloatMatrix& u, const FloatMatrix& v)
545 {
546 if (packed ())
547 unpack ();
548
549 FloatMatrix& l = m_L;
550 FloatMatrix& r = m_a_fact;
551
552 F77_INT m = to_f77_int (l.rows ());
553 F77_INT n = to_f77_int (r.columns ());
554 F77_INT k = to_f77_int (l.columns ());
555
556 F77_INT u_nr = to_f77_int (u.rows ());
557 F77_INT u_nc = to_f77_int (u.columns ());
558
559 F77_INT v_nr = to_f77_int (v.rows ());
560 F77_INT v_nc = to_f77_int (v.columns ());
561
562 if (u_nr != m || v_nr != n || u_nc != v_nc)
563 (*current_liboctave_error_handler) ("luupdate: dimensions mismatch");
564
565 OCTAVE_LOCAL_BUFFER (float, w, m);
566 for (F77_INT i = 0; i < m; i++) m_ipvt(i) += 1; // increment
567 for (volatile F77_INT i = 0; i < u_nc; i++)
568 {
569 FloatColumnVector utmp = u.column (i);
570 FloatColumnVector vtmp = v.column (i);
571 F77_XFCN (slup1up, SLUP1UP, (m, n, l.fortran_vec (),
572 m, r.fortran_vec (), k,
573 m_ipvt.fortran_vec (),
574 utmp.data (), vtmp.data (), w));
575 }
576 for (F77_INT i = 0; i < m; i++) m_ipvt(i) -= 1; // decrement
577 }
578 578
579 #endif 579 #endif
580 580
581 template <> 581 template <>
582 OCTAVE_API 582 OCTAVE_API
583 lu<ComplexMatrix>::lu (const ComplexMatrix& a) 583 lu<ComplexMatrix>::lu (const ComplexMatrix& a)
584 { 584 {
585 F77_INT a_nr = to_f77_int (a.rows ()); 585 F77_INT a_nr = to_f77_int (a.rows ());
586 F77_INT a_nc = to_f77_int (a.columns ()); 586 F77_INT a_nc = to_f77_int (a.columns ());
587 F77_INT mn = (a_nr < a_nc ? a_nr : a_nc); 587 F77_INT mn = (a_nr < a_nc ? a_nr : a_nc);
588 588
589 m_ipvt.resize (dim_vector (mn, 1)); 589 m_ipvt.resize (dim_vector (mn, 1));
590 F77_INT *pipvt = m_ipvt.fortran_vec (); 590 F77_INT *pipvt = m_ipvt.fortran_vec ();
591 591
592 m_a_fact = a; 592 m_a_fact = a;
593 Complex *tmp_data = m_a_fact.fortran_vec (); 593 Complex *tmp_data = m_a_fact.fortran_vec ();
594 594
595 F77_INT info = 0; 595 F77_INT info = 0;
596 596
597 F77_XFCN (zgetrf, ZGETRF, (a_nr, a_nc, F77_DBLE_CMPLX_ARG (tmp_data), 597 F77_XFCN (zgetrf, ZGETRF, (a_nr, a_nc, F77_DBLE_CMPLX_ARG (tmp_data),
598 a_nr, pipvt, info)); 598 a_nr, pipvt, info));
599 599
600 for (F77_INT i = 0; i < mn; i++) 600 for (F77_INT i = 0; i < mn; i++)
601 pipvt[i] -= 1; 601 pipvt[i] -= 1;
602 } 602 }
603 603
604 #if defined (HAVE_QRUPDATE_LUU) 604 #if defined (HAVE_QRUPDATE_LUU)
605 605
606 template <> 606 template <>
607 OCTAVE_API void 607 OCTAVE_API void
608 lu<ComplexMatrix>::update (const ComplexColumnVector& u, 608 lu<ComplexMatrix>::update (const ComplexColumnVector& u,
609 const ComplexColumnVector& v) 609 const ComplexColumnVector& v)
610 { 610 {
611 if (packed ()) 611 if (packed ())
612 unpack (); 612 unpack ();
613 613
614 ComplexMatrix& l = m_L; 614 ComplexMatrix& l = m_L;
615 ComplexMatrix& r = m_a_fact; 615 ComplexMatrix& r = m_a_fact;
616 616
617 F77_INT m = to_f77_int (l.rows ()); 617 F77_INT m = to_f77_int (l.rows ());
618 F77_INT n = to_f77_int (r.columns ()); 618 F77_INT n = to_f77_int (r.columns ());
619 F77_INT k = to_f77_int (l.columns ()); 619 F77_INT k = to_f77_int (l.columns ());
620 620
621 F77_INT u_nel = to_f77_int (u.numel ()); 621 F77_INT u_nel = to_f77_int (u.numel ());
622 F77_INT v_nel = to_f77_int (v.numel ()); 622 F77_INT v_nel = to_f77_int (v.numel ());
623 623
624 if (u_nel != m || v_nel != n) 624 if (u_nel != m || v_nel != n)
625 (*current_liboctave_error_handler) ("luupdate: dimensions mismatch"); 625 (*current_liboctave_error_handler) ("luupdate: dimensions mismatch");
626 626
627 ComplexColumnVector utmp = u; 627 ComplexColumnVector utmp = u;
628 ComplexColumnVector vtmp = v; 628 ComplexColumnVector vtmp = v;
629 F77_XFCN (zlu1up, ZLU1UP, (m, n, F77_DBLE_CMPLX_ARG (l.fortran_vec ()), m, 629 F77_XFCN (zlu1up, ZLU1UP, (m, n, F77_DBLE_CMPLX_ARG (l.fortran_vec ()), m,
630 F77_DBLE_CMPLX_ARG (r.fortran_vec ()), k, 630 F77_DBLE_CMPLX_ARG (r.fortran_vec ()), k,
631 F77_DBLE_CMPLX_ARG (utmp.fortran_vec ()),
632 F77_DBLE_CMPLX_ARG (vtmp.fortran_vec ())));
633 }
634
635 template <>
636 OCTAVE_API void
637 lu<ComplexMatrix>::update (const ComplexMatrix& u, const ComplexMatrix& v)
638 {
639 if (packed ())
640 unpack ();
641
642 ComplexMatrix& l = m_L;
643 ComplexMatrix& r = m_a_fact;
644
645 F77_INT m = to_f77_int (l.rows ());
646 F77_INT n = to_f77_int (r.columns ());
647 F77_INT k = to_f77_int (l.columns ());
648
649 F77_INT u_nr = to_f77_int (u.rows ());
650 F77_INT u_nc = to_f77_int (u.columns ());
651
652 F77_INT v_nr = to_f77_int (v.rows ());
653 F77_INT v_nc = to_f77_int (v.columns ());
654
655 if (u_nr != m || v_nr != n || u_nc != v_nc)
656 (*current_liboctave_error_handler) ("luupdate: dimensions mismatch");
657
658 for (volatile F77_INT i = 0; i < u_nc; i++)
659 {
660 ComplexColumnVector utmp = u.column (i);
661 ComplexColumnVector vtmp = v.column (i);
662 F77_XFCN (zlu1up, ZLU1UP, (m, n,
663 F77_DBLE_CMPLX_ARG (l.fortran_vec ()),
664 m, F77_DBLE_CMPLX_ARG (r.fortran_vec ()),
665 k,
631 F77_DBLE_CMPLX_ARG (utmp.fortran_vec ()), 666 F77_DBLE_CMPLX_ARG (utmp.fortran_vec ()),
632 F77_DBLE_CMPLX_ARG (vtmp.fortran_vec ()))); 667 F77_DBLE_CMPLX_ARG (vtmp.fortran_vec ())));
633 } 668 }
634 669 }
635 template <> 670
636 OCTAVE_API void 671 template <>
637 lu<ComplexMatrix>::update (const ComplexMatrix& u, const ComplexMatrix& v) 672 OCTAVE_API void
638 { 673 lu<ComplexMatrix>::update_piv (const ComplexColumnVector& u,
639 if (packed ()) 674 const ComplexColumnVector& v)
640 unpack (); 675 {
641 676 if (packed ())
642 ComplexMatrix& l = m_L; 677 unpack ();
643 ComplexMatrix& r = m_a_fact; 678
644 679 ComplexMatrix& l = m_L;
645 F77_INT m = to_f77_int (l.rows ()); 680 ComplexMatrix& r = m_a_fact;
646 F77_INT n = to_f77_int (r.columns ()); 681
647 F77_INT k = to_f77_int (l.columns ()); 682 F77_INT m = to_f77_int (l.rows ());
648 683 F77_INT n = to_f77_int (r.columns ());
649 F77_INT u_nr = to_f77_int (u.rows ()); 684 F77_INT k = to_f77_int (l.columns ());
650 F77_INT u_nc = to_f77_int (u.columns ()); 685
651 686 F77_INT u_nel = to_f77_int (u.numel ());
652 F77_INT v_nr = to_f77_int (v.rows ()); 687 F77_INT v_nel = to_f77_int (v.numel ());
653 F77_INT v_nc = to_f77_int (v.columns ()); 688
654 689 if (u_nel != m || v_nel != n)
655 if (u_nr != m || v_nr != n || u_nc != v_nc) 690 (*current_liboctave_error_handler) ("luupdate: dimensions mismatch");
656 (*current_liboctave_error_handler) ("luupdate: dimensions mismatch"); 691
657 692 ComplexColumnVector utmp = u;
658 for (volatile F77_INT i = 0; i < u_nc; i++) 693 ComplexColumnVector vtmp = v;
659 { 694 OCTAVE_LOCAL_BUFFER (Complex, w, m);
660 ComplexColumnVector utmp = u.column (i); 695 for (F77_INT i = 0; i < m; i++) m_ipvt(i) += 1; // increment
661 ComplexColumnVector vtmp = v.column (i); 696 F77_XFCN (zlup1up, ZLUP1UP, (m, n, F77_DBLE_CMPLX_ARG (l.fortran_vec ()),
662 F77_XFCN (zlu1up, ZLU1UP, (m, n, 697 m, F77_DBLE_CMPLX_ARG (r.fortran_vec ()), k,
663 F77_DBLE_CMPLX_ARG (l.fortran_vec ()), 698 m_ipvt.fortran_vec (),
664 m, F77_DBLE_CMPLX_ARG (r.fortran_vec ()), 699 F77_CONST_DBLE_CMPLX_ARG (utmp.data ()),
665 k, 700 F77_CONST_DBLE_CMPLX_ARG (vtmp.data ()),
666 F77_DBLE_CMPLX_ARG (utmp.fortran_vec ()), 701 F77_DBLE_CMPLX_ARG (w)));
667 F77_DBLE_CMPLX_ARG (vtmp.fortran_vec ()))); 702 for (F77_INT i = 0; i < m; i++) m_ipvt(i) -= 1; // decrement
668 } 703 }
669 } 704
670 705 template <>
671 template <> 706 OCTAVE_API void
672 OCTAVE_API void 707 lu<ComplexMatrix>::update_piv (const ComplexMatrix& u,
673 lu<ComplexMatrix>::update_piv (const ComplexColumnVector& u, 708 const ComplexMatrix& v)
674 const ComplexColumnVector& v) 709 {
675 { 710 if (packed ())
676 if (packed ()) 711 unpack ();
677 unpack (); 712
678 713 ComplexMatrix& l = m_L;
679 ComplexMatrix& l = m_L; 714 ComplexMatrix& r = m_a_fact;
680 ComplexMatrix& r = m_a_fact; 715
681 716 F77_INT m = to_f77_int (l.rows ());
682 F77_INT m = to_f77_int (l.rows ()); 717 F77_INT n = to_f77_int (r.columns ());
683 F77_INT n = to_f77_int (r.columns ()); 718 F77_INT k = to_f77_int (l.columns ());
684 F77_INT k = to_f77_int (l.columns ()); 719
685 720 F77_INT u_nr = to_f77_int (u.rows ());
686 F77_INT u_nel = to_f77_int (u.numel ()); 721 F77_INT u_nc = to_f77_int (u.columns ());
687 F77_INT v_nel = to_f77_int (v.numel ()); 722
688 723 F77_INT v_nr = to_f77_int (v.rows ());
689 if (u_nel != m || v_nel != n) 724 F77_INT v_nc = to_f77_int (v.columns ());
690 (*current_liboctave_error_handler) ("luupdate: dimensions mismatch"); 725
691 726 if (u_nr != m || v_nr != n || u_nc != v_nc)
692 ComplexColumnVector utmp = u; 727 (*current_liboctave_error_handler) ("luupdate: dimensions mismatch");
693 ComplexColumnVector vtmp = v; 728
694 OCTAVE_LOCAL_BUFFER (Complex, w, m); 729 OCTAVE_LOCAL_BUFFER (Complex, w, m);
695 for (F77_INT i = 0; i < m; i++) m_ipvt(i) += 1; // increment 730 for (F77_INT i = 0; i < m; i++) m_ipvt(i) += 1; // increment
696 F77_XFCN (zlup1up, ZLUP1UP, (m, n, F77_DBLE_CMPLX_ARG (l.fortran_vec ()), 731 for (volatile F77_INT i = 0; i < u_nc; i++)
697 m, F77_DBLE_CMPLX_ARG (r.fortran_vec ()), k, 732 {
698 m_ipvt.fortran_vec (), 733 ComplexColumnVector utmp = u.column (i);
734 ComplexColumnVector vtmp = v.column (i);
735 F77_XFCN (zlup1up, ZLUP1UP, (m, n,
736 F77_DBLE_CMPLX_ARG (l.fortran_vec ()),
737 m,
738 F77_DBLE_CMPLX_ARG (r.fortran_vec ()),
739 k, m_ipvt.fortran_vec (),
699 F77_CONST_DBLE_CMPLX_ARG (utmp.data ()), 740 F77_CONST_DBLE_CMPLX_ARG (utmp.data ()),
700 F77_CONST_DBLE_CMPLX_ARG (vtmp.data ()), 741 F77_CONST_DBLE_CMPLX_ARG (vtmp.data ()),
701 F77_DBLE_CMPLX_ARG (w))); 742 F77_DBLE_CMPLX_ARG (w)));
702 for (F77_INT i = 0; i < m; i++) m_ipvt(i) -= 1; // decrement 743 }
703 } 744 for (F77_INT i = 0; i < m; i++) m_ipvt(i) -= 1; // decrement
704 745 }
705 template <>
706 OCTAVE_API void
707 lu<ComplexMatrix>::update_piv (const ComplexMatrix& u,
708 const ComplexMatrix& v)
709 {
710 if (packed ())
711 unpack ();
712
713 ComplexMatrix& l = m_L;
714 ComplexMatrix& r = m_a_fact;
715
716 F77_INT m = to_f77_int (l.rows ());
717 F77_INT n = to_f77_int (r.columns ());
718 F77_INT k = to_f77_int (l.columns ());
719
720 F77_INT u_nr = to_f77_int (u.rows ());
721 F77_INT u_nc = to_f77_int (u.columns ());
722
723 F77_INT v_nr = to_f77_int (v.rows ());
724 F77_INT v_nc = to_f77_int (v.columns ());
725
726 if (u_nr != m || v_nr != n || u_nc != v_nc)
727 (*current_liboctave_error_handler) ("luupdate: dimensions mismatch");
728
729 OCTAVE_LOCAL_BUFFER (Complex, w, m);
730 for (F77_INT i = 0; i < m; i++) m_ipvt(i) += 1; // increment
731 for (volatile F77_INT i = 0; i < u_nc; i++)
732 {
733 ComplexColumnVector utmp = u.column (i);
734 ComplexColumnVector vtmp = v.column (i);
735 F77_XFCN (zlup1up, ZLUP1UP, (m, n,
736 F77_DBLE_CMPLX_ARG (l.fortran_vec ()),
737 m,
738 F77_DBLE_CMPLX_ARG (r.fortran_vec ()),
739 k, m_ipvt.fortran_vec (),
740 F77_CONST_DBLE_CMPLX_ARG (utmp.data ()),
741 F77_CONST_DBLE_CMPLX_ARG (vtmp.data ()),
742 F77_DBLE_CMPLX_ARG (w)));
743 }
744 for (F77_INT i = 0; i < m; i++) m_ipvt(i) -= 1; // decrement
745 }
746 746
747 #endif 747 #endif
748 748
749 template <> 749 template <>
750 OCTAVE_API 750 OCTAVE_API
751 lu<FloatComplexMatrix>::lu (const FloatComplexMatrix& a) 751 lu<FloatComplexMatrix>::lu (const FloatComplexMatrix& a)
752 { 752 {
753 F77_INT a_nr = to_f77_int (a.rows ()); 753 F77_INT a_nr = to_f77_int (a.rows ());
754 F77_INT a_nc = to_f77_int (a.columns ()); 754 F77_INT a_nc = to_f77_int (a.columns ());
755 F77_INT mn = (a_nr < a_nc ? a_nr : a_nc); 755 F77_INT mn = (a_nr < a_nc ? a_nr : a_nc);
756 756
757 m_ipvt.resize (dim_vector (mn, 1)); 757 m_ipvt.resize (dim_vector (mn, 1));
758 F77_INT *pipvt = m_ipvt.fortran_vec (); 758 F77_INT *pipvt = m_ipvt.fortran_vec ();
759 759
760 m_a_fact = a; 760 m_a_fact = a;
761 FloatComplex *tmp_data = m_a_fact.fortran_vec (); 761 FloatComplex *tmp_data = m_a_fact.fortran_vec ();
762 762
763 F77_INT info = 0; 763 F77_INT info = 0;
764 764
765 F77_XFCN (cgetrf, CGETRF, (a_nr, a_nc, F77_CMPLX_ARG (tmp_data), a_nr, 765 F77_XFCN (cgetrf, CGETRF, (a_nr, a_nc, F77_CMPLX_ARG (tmp_data), a_nr,
766 pipvt, info)); 766 pipvt, info));
767 767
768 for (F77_INT i = 0; i < mn; i++) 768 for (F77_INT i = 0; i < mn; i++)
769 pipvt[i] -= 1; 769 pipvt[i] -= 1;
770 } 770 }
771 771
772 #if defined (HAVE_QRUPDATE_LUU) 772 #if defined (HAVE_QRUPDATE_LUU)
773 773
774 template <> 774 template <>
775 OCTAVE_API void 775 OCTAVE_API void
776 lu<FloatComplexMatrix>::update (const FloatComplexColumnVector& u, 776 lu<FloatComplexMatrix>::update (const FloatComplexColumnVector& u,
777 const FloatComplexColumnVector& v) 777 const FloatComplexColumnVector& v)
778 { 778 {
779 if (packed ()) 779 if (packed ())
780 unpack (); 780 unpack ();
781 781
782 FloatComplexMatrix& l = m_L; 782 FloatComplexMatrix& l = m_L;
783 FloatComplexMatrix& r = m_a_fact; 783 FloatComplexMatrix& r = m_a_fact;
784 784
785 F77_INT m = to_f77_int (l.rows ()); 785 F77_INT m = to_f77_int (l.rows ());
786 F77_INT n = to_f77_int (r.columns ()); 786 F77_INT n = to_f77_int (r.columns ());
787 F77_INT k = to_f77_int (l.columns ()); 787 F77_INT k = to_f77_int (l.columns ());
788 788
789 F77_INT u_nel = to_f77_int (u.numel ()); 789 F77_INT u_nel = to_f77_int (u.numel ());
790 F77_INT v_nel = to_f77_int (v.numel ()); 790 F77_INT v_nel = to_f77_int (v.numel ());
791 791
792 if (u_nel != m || v_nel != n) 792 if (u_nel != m || v_nel != n)
793 (*current_liboctave_error_handler) ("luupdate: dimensions mismatch"); 793 (*current_liboctave_error_handler) ("luupdate: dimensions mismatch");
794 794
795 FloatComplexColumnVector utmp = u; 795 FloatComplexColumnVector utmp = u;
796 FloatComplexColumnVector vtmp = v; 796 FloatComplexColumnVector vtmp = v;
797 F77_XFCN (clu1up, CLU1UP, (m, n, F77_CMPLX_ARG (l.fortran_vec ()), m, 797 F77_XFCN (clu1up, CLU1UP, (m, n, F77_CMPLX_ARG (l.fortran_vec ()), m,
798 F77_CMPLX_ARG (r.fortran_vec ()), k, 798 F77_CMPLX_ARG (r.fortran_vec ()), k,
799 F77_CMPLX_ARG (utmp.fortran_vec ()),
800 F77_CMPLX_ARG (vtmp.fortran_vec ())));
801 }
802
803 template <>
804 OCTAVE_API void
805 lu<FloatComplexMatrix>::update (const FloatComplexMatrix& u,
806 const FloatComplexMatrix& v)
807 {
808 if (packed ())
809 unpack ();
810
811 FloatComplexMatrix& l = m_L;
812 FloatComplexMatrix& r = m_a_fact;
813
814 F77_INT m = to_f77_int (l.rows ());
815 F77_INT n = to_f77_int (r.columns ());
816 F77_INT k = to_f77_int (l.columns ());
817
818 F77_INT u_nr = to_f77_int (u.rows ());
819 F77_INT u_nc = to_f77_int (u.columns ());
820
821 F77_INT v_nr = to_f77_int (v.rows ());
822 F77_INT v_nc = to_f77_int (v.columns ());
823
824 if (u_nr != m || v_nr != n || u_nc != v_nc)
825 (*current_liboctave_error_handler) ("luupdate: dimensions mismatch");
826
827 for (volatile F77_INT i = 0; i < u_nc; i++)
828 {
829 FloatComplexColumnVector utmp = u.column (i);
830 FloatComplexColumnVector vtmp = v.column (i);
831 F77_XFCN (clu1up, CLU1UP, (m, n, F77_CMPLX_ARG (l.fortran_vec ()),
832 m, F77_CMPLX_ARG (r.fortran_vec ()), k,
799 F77_CMPLX_ARG (utmp.fortran_vec ()), 833 F77_CMPLX_ARG (utmp.fortran_vec ()),
800 F77_CMPLX_ARG (vtmp.fortran_vec ()))); 834 F77_CMPLX_ARG (vtmp.fortran_vec ())));
801 } 835 }
802 836 }
803 template <> 837
804 OCTAVE_API void 838 template <>
805 lu<FloatComplexMatrix>::update (const FloatComplexMatrix& u, 839 OCTAVE_API void
840 lu<FloatComplexMatrix>::update_piv (const FloatComplexColumnVector& u,
841 const FloatComplexColumnVector& v)
842 {
843 if (packed ())
844 unpack ();
845
846 FloatComplexMatrix& l = m_L;
847 FloatComplexMatrix& r = m_a_fact;
848
849 F77_INT m = to_f77_int (l.rows ());
850 F77_INT n = to_f77_int (r.columns ());
851 F77_INT k = to_f77_int (l.columns ());
852
853 F77_INT u_nel = to_f77_int (u.numel ());
854 F77_INT v_nel = to_f77_int (v.numel ());
855
856 if (u_nel != m || v_nel != n)
857 (*current_liboctave_error_handler) ("luupdate: dimensions mismatch");
858
859 FloatComplexColumnVector utmp = u;
860 FloatComplexColumnVector vtmp = v;
861 OCTAVE_LOCAL_BUFFER (FloatComplex, w, m);
862 for (F77_INT i = 0; i < m; i++) m_ipvt(i) += 1; // increment
863 F77_XFCN (clup1up, CLUP1UP, (m, n, F77_CMPLX_ARG (l.fortran_vec ()),
864 m, F77_CMPLX_ARG (r.fortran_vec ()), k,
865 m_ipvt.fortran_vec (),
866 F77_CONST_CMPLX_ARG (utmp.data ()),
867 F77_CONST_CMPLX_ARG (vtmp.data ()),
868 F77_CMPLX_ARG (w)));
869 for (F77_INT i = 0; i < m; i++) m_ipvt(i) -= 1; // decrement
870 }
871
872 template <>
873 OCTAVE_API void
874 lu<FloatComplexMatrix>::update_piv (const FloatComplexMatrix& u,
806 const FloatComplexMatrix& v) 875 const FloatComplexMatrix& v)
807 { 876 {
808 if (packed ()) 877 if (packed ())
809 unpack (); 878 unpack ();
810 879
811 FloatComplexMatrix& l = m_L; 880 FloatComplexMatrix& l = m_L;
812 FloatComplexMatrix& r = m_a_fact; 881 FloatComplexMatrix& r = m_a_fact;
813 882
814 F77_INT m = to_f77_int (l.rows ()); 883 F77_INT m = to_f77_int (l.rows ());
815 F77_INT n = to_f77_int (r.columns ()); 884 F77_INT n = to_f77_int (r.columns ());
816 F77_INT k = to_f77_int (l.columns ()); 885 F77_INT k = to_f77_int (l.columns ());
817 886
818 F77_INT u_nr = to_f77_int (u.rows ()); 887 F77_INT u_nr = to_f77_int (u.rows ());
819 F77_INT u_nc = to_f77_int (u.columns ()); 888 F77_INT u_nc = to_f77_int (u.columns ());
820 889
821 F77_INT v_nr = to_f77_int (v.rows ()); 890 F77_INT v_nr = to_f77_int (v.rows ());
822 F77_INT v_nc = to_f77_int (v.columns ()); 891 F77_INT v_nc = to_f77_int (v.columns ());
823 892
824 if (u_nr != m || v_nr != n || u_nc != v_nc) 893 if (u_nr != m || v_nr != n || u_nc != v_nc)
825 (*current_liboctave_error_handler) ("luupdate: dimensions mismatch"); 894 (*current_liboctave_error_handler) ("luupdate: dimensions mismatch");
826 895
827 for (volatile F77_INT i = 0; i < u_nc; i++) 896 OCTAVE_LOCAL_BUFFER (FloatComplex, w, m);
828 { 897 for (F77_INT i = 0; i < m; i++) m_ipvt(i) += 1; // increment
829 FloatComplexColumnVector utmp = u.column (i); 898 for (volatile F77_INT i = 0; i < u_nc; i++)
830 FloatComplexColumnVector vtmp = v.column (i); 899 {
831 F77_XFCN (clu1up, CLU1UP, (m, n, F77_CMPLX_ARG (l.fortran_vec ()), 900 FloatComplexColumnVector utmp = u.column (i);
832 m, F77_CMPLX_ARG (r.fortran_vec ()), k, 901 FloatComplexColumnVector vtmp = v.column (i);
833 F77_CMPLX_ARG (utmp.fortran_vec ()),
834 F77_CMPLX_ARG (vtmp.fortran_vec ())));
835 }
836 }
837
838 template <>
839 OCTAVE_API void
840 lu<FloatComplexMatrix>::update_piv (const FloatComplexColumnVector& u,
841 const FloatComplexColumnVector& v)
842 {
843 if (packed ())
844 unpack ();
845
846 FloatComplexMatrix& l = m_L;
847 FloatComplexMatrix& r = m_a_fact;
848
849 F77_INT m = to_f77_int (l.rows ());
850 F77_INT n = to_f77_int (r.columns ());
851 F77_INT k = to_f77_int (l.columns ());
852
853 F77_INT u_nel = to_f77_int (u.numel ());
854 F77_INT v_nel = to_f77_int (v.numel ());
855
856 if (u_nel != m || v_nel != n)
857 (*current_liboctave_error_handler) ("luupdate: dimensions mismatch");
858
859 FloatComplexColumnVector utmp = u;
860 FloatComplexColumnVector vtmp = v;
861 OCTAVE_LOCAL_BUFFER (FloatComplex, w, m);
862 for (F77_INT i = 0; i < m; i++) m_ipvt(i) += 1; // increment
863 F77_XFCN (clup1up, CLUP1UP, (m, n, F77_CMPLX_ARG (l.fortran_vec ()), 902 F77_XFCN (clup1up, CLUP1UP, (m, n, F77_CMPLX_ARG (l.fortran_vec ()),
864 m, F77_CMPLX_ARG (r.fortran_vec ()), k, 903 m, F77_CMPLX_ARG (r.fortran_vec ()), k,
865 m_ipvt.fortran_vec (), 904 m_ipvt.fortran_vec (),
866 F77_CONST_CMPLX_ARG (utmp.data ()), 905 F77_CONST_CMPLX_ARG (utmp.data ()),
867 F77_CONST_CMPLX_ARG (vtmp.data ()), 906 F77_CONST_CMPLX_ARG (vtmp.data ()),
868 F77_CMPLX_ARG (w))); 907 F77_CMPLX_ARG (w)));
869 for (F77_INT i = 0; i < m; i++) m_ipvt(i) -= 1; // decrement 908 }
870 } 909 for (F77_INT i = 0; i < m; i++) m_ipvt(i) -= 1; // decrement
871 910 }
872 template <>
873 OCTAVE_API void
874 lu<FloatComplexMatrix>::update_piv (const FloatComplexMatrix& u,
875 const FloatComplexMatrix& v)
876 {
877 if (packed ())
878 unpack ();
879
880 FloatComplexMatrix& l = m_L;
881 FloatComplexMatrix& r = m_a_fact;
882
883 F77_INT m = to_f77_int (l.rows ());
884 F77_INT n = to_f77_int (r.columns ());
885 F77_INT k = to_f77_int (l.columns ());
886
887 F77_INT u_nr = to_f77_int (u.rows ());
888 F77_INT u_nc = to_f77_int (u.columns ());
889
890 F77_INT v_nr = to_f77_int (v.rows ());
891 F77_INT v_nc = to_f77_int (v.columns ());
892
893 if (u_nr != m || v_nr != n || u_nc != v_nc)
894 (*current_liboctave_error_handler) ("luupdate: dimensions mismatch");
895
896 OCTAVE_LOCAL_BUFFER (FloatComplex, w, m);
897 for (F77_INT i = 0; i < m; i++) m_ipvt(i) += 1; // increment
898 for (volatile F77_INT i = 0; i < u_nc; i++)
899 {
900 FloatComplexColumnVector utmp = u.column (i);
901 FloatComplexColumnVector vtmp = v.column (i);
902 F77_XFCN (clup1up, CLUP1UP, (m, n, F77_CMPLX_ARG (l.fortran_vec ()),
903 m, F77_CMPLX_ARG (r.fortran_vec ()), k,
904 m_ipvt.fortran_vec (),
905 F77_CONST_CMPLX_ARG (utmp.data ()),
906 F77_CONST_CMPLX_ARG (vtmp.data ()),
907 F77_CMPLX_ARG (w)));
908 }
909 for (F77_INT i = 0; i < m; i++) m_ipvt(i) -= 1; // decrement
910 }
911 911
912 #endif 912 #endif
913 913
914 // Instantiations we need. 914 // Instantiations we need.
915 915
916 template class lu<Matrix>; 916 template class lu<Matrix>;
917 917
918 template class lu<FloatMatrix>; 918 template class lu<FloatMatrix>;
919 919
920 template class lu<ComplexMatrix>; 920 template class lu<ComplexMatrix>;
921 921
922 template class lu<FloatComplexMatrix>; 922 template class lu<FloatComplexMatrix>;
923 923
924 OCTAVE_END_NAMESPACE(math) 924 OCTAVE_END_NAMESPACE(math)
925 OCTAVE_END_NAMESPACE(octave) 925 OCTAVE_END_NAMESPACE(octave)