GNU Octave  9.1.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
lu.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 //
3 // Copyright (C) 1996-2024 The Octave Project Developers
4 //
5 // See the file COPYRIGHT.md in the top-level directory of this
6 // distribution or <https://octave.org/copyright/>.
7 //
8 // This file is part of Octave.
9 //
10 // Octave is free software: you can redistribute it and/or modify it
11 // under the terms of the GNU General Public License as published by
12 // the Free Software Foundation, either version 3 of the License, or
13 // (at your option) any later version.
14 //
15 // Octave is distributed in the hope that it will be useful, but
16 // WITHOUT ANY WARRANTY; without even the implied warranty of
17 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18 // GNU General Public License for more details.
19 //
20 // You should have received a copy of the GNU General Public License
21 // along with Octave; see the file COPYING. If not, see
22 // <https://www.gnu.org/licenses/>.
23 //
24 ////////////////////////////////////////////////////////////////////////
25 
26 #if defined (HAVE_CONFIG_H)
27 # include "config.h"
28 #endif
29 
30 #include <algorithm>
31 
32 #include "CColVector.h"
33 #include "CMatrix.h"
34 #include "PermMatrix.h"
35 #include "dColVector.h"
36 #include "dMatrix.h"
37 #include "fCColVector.h"
38 #include "fCMatrix.h"
39 #include "fColVector.h"
40 #include "fMatrix.h"
41 #include "lo-error.h"
42 #include "lo-lapack-proto.h"
43 #include "lo-qrupdate-proto.h"
44 #include "lu.h"
45 #include "oct-locbuf.h"
46 
47 OCTAVE_BEGIN_NAMESPACE(octave)
48 
49 OCTAVE_BEGIN_NAMESPACE(math)
50 
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
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
55 // Fluupdate, it could be given values that are out of range. We
56 // should ensure that the values are within range here.
57 
58 template <typename T>
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 ())
61 {
62  if (l.columns () != u.rows ())
63  (*current_liboctave_error_handler) ("lu: dimension mismatch");
64 }
65 
66 template <typename T>
67 bool
68 lu<T>::packed () const
69 {
70  return m_L.dims () == dim_vector ();
71 }
72 
73 template <typename T>
74 void
75 lu<T>::unpack ()
76 {
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 () 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++)
107  {
108  if (i < a_nc)
109  l.xelem (i, i) = 1.0;
110 
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 () 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 () 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 () 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 () const
190 {
191  return PermMatrix (getp (), false);
192 }
193 
194 template <typename T>
195 ColumnVector
196 lu<T>::P_vec () 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 () 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 
230 #if ! defined (HAVE_QRUPDATE_LUU)
231 
232 template <typename T>
233 void
234 lu<T>::update (const VT&, const VT&)
235 {
236  (*current_liboctave_error_handler)
237  ("luupdate: support for qrupdate with LU updates "
238  "was unavailable or disabled when liboctave was built");
239 }
240 
241 template <typename T>
242 void
243 lu<T>::update (const T&, const T&)
244 {
245  (*current_liboctave_error_handler)
246  ("luupdate: support for qrupdate with LU updates "
247  "was unavailable or disabled when liboctave was built");
248 }
249 
250 template <typename T>
251 void
252 lu<T>::update_piv (const VT&, const VT&)
253 {
254  (*current_liboctave_error_handler)
255  ("luupdate: support for qrupdate with LU updates "
256  "was unavailable or disabled when liboctave was built");
257 }
258 
259 template <typename T>
260 void
261 lu<T>::update_piv (const T&, const T&)
262 {
263  (*current_liboctave_error_handler)
264  ("luupdate: support for qrupdate with LU updates "
265  "was unavailable or disabled when liboctave was built");
266 }
267 
268 #endif
269 
270 // Specializations.
271 
272 template <>
273 OCTAVE_API
274 lu<Matrix>::lu (const Matrix& a)
275 {
276  F77_INT a_nr = to_f77_int (a.rows ());
277  F77_INT a_nc = to_f77_int (a.columns ());
278  F77_INT mn = (a_nr < a_nc ? a_nr : a_nc);
279 
280  m_ipvt.resize (dim_vector (mn, 1));
281  F77_INT *pipvt = m_ipvt.fortran_vec ();
282 
283  m_a_fact = a;
284  double *tmp_data = m_a_fact.fortran_vec ();
285 
286  F77_INT info = 0;
287 
288  F77_XFCN (dgetrf, DGETRF, (a_nr, a_nc, tmp_data, a_nr, pipvt, info));
289 
290  for (F77_INT i = 0; i < mn; i++)
291  pipvt[i] -= 1;
292 }
293 
294 #if defined (HAVE_QRUPDATE_LUU)
295 
296 template <>
297 OCTAVE_API void
298 lu<Matrix>::update (const ColumnVector& u, const ColumnVector& v)
299 {
300  if (packed ())
301  unpack ();
302 
303  Matrix& l = m_L;
304  Matrix& r = m_a_fact;
305 
306  F77_INT m = to_f77_int (l.rows ());
307  F77_INT n = to_f77_int (r.columns ());
308  F77_INT k = to_f77_int (l.columns ());
309 
310  F77_INT u_nel = to_f77_int (u.numel ());
311  F77_INT v_nel = to_f77_int (v.numel ());
312 
313  if (u_nel != m || v_nel != n)
314  (*current_liboctave_error_handler) ("luupdate: dimensions mismatch");
315 
316  ColumnVector utmp = u;
317  ColumnVector vtmp = v;
318  F77_XFCN (dlu1up, DLU1UP, (m, n, l.fortran_vec (), m, r.fortran_vec (),
319  k, utmp.fortran_vec (), vtmp.fortran_vec ()));
320 }
321 
322 template <>
323 OCTAVE_API void
324 lu<Matrix>::update (const Matrix& u, const Matrix& v)
325 {
326  if (packed ())
327  unpack ();
328 
329  Matrix& l = m_L;
330  Matrix& r = m_a_fact;
331 
332  F77_INT m = to_f77_int (l.rows ());
333  F77_INT n = to_f77_int (r.columns ());
334  F77_INT k = to_f77_int (l.columns ());
335 
336  F77_INT u_nr = to_f77_int (u.rows ());
337  F77_INT u_nc = to_f77_int (u.columns ());
338 
339  F77_INT v_nr = to_f77_int (v.rows ());
340  F77_INT v_nc = to_f77_int (v.columns ());
341 
342  if (u_nr != m || v_nr != n || u_nc != v_nc)
343  (*current_liboctave_error_handler) ("luupdate: dimensions mismatch");
344 
345  for (volatile F77_INT i = 0; i < u_nc; i++)
346  {
347  ColumnVector utmp = u.column (i);
348  ColumnVector vtmp = v.column (i);
349  F77_XFCN (dlu1up, DLU1UP, (m, n, l.fortran_vec (),
350  m, r.fortran_vec (), k,
351  utmp.fortran_vec (), vtmp.fortran_vec ()));
352  }
353 }
354 
355 template <>
356 OCTAVE_API void
357 lu<Matrix>::update_piv (const ColumnVector& u, const ColumnVector& v)
358 {
359  if (packed ())
360  unpack ();
361 
362  Matrix& l = m_L;
363  Matrix& r = m_a_fact;
364 
365  F77_INT m = to_f77_int (l.rows ());
366  F77_INT n = to_f77_int (r.columns ());
367  F77_INT k = to_f77_int (l.columns ());
368 
369  F77_INT u_nel = to_f77_int (u.numel ());
370  F77_INT v_nel = to_f77_int (v.numel ());
371 
372  if (u_nel != m || v_nel != n)
373  (*current_liboctave_error_handler) ("luupdate: dimensions mismatch");
374 
375  ColumnVector utmp = u;
376  ColumnVector vtmp = v;
377  OCTAVE_LOCAL_BUFFER (double, w, m);
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);
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 
423 #endif
424 
425 template <>
426 OCTAVE_API
427 lu<FloatMatrix>::lu (const FloatMatrix& a)
428 {
429  F77_INT a_nr = to_f77_int (a.rows ());
430  F77_INT a_nc = to_f77_int (a.columns ());
431  F77_INT mn = (a_nr < a_nc ? a_nr : a_nc);
432 
433  m_ipvt.resize (dim_vector (mn, 1));
434  F77_INT *pipvt = m_ipvt.fortran_vec ();
435 
436  m_a_fact = a;
437  float *tmp_data = m_a_fact.fortran_vec ();
438 
439  F77_INT info = 0;
440 
441  F77_XFCN (sgetrf, SGETRF, (a_nr, a_nc, tmp_data, a_nr, pipvt, info));
442 
443  for (F77_INT i = 0; i < mn; i++)
444  pipvt[i] -= 1;
445 }
446 
447 #if defined (HAVE_QRUPDATE_LUU)
448 
449 template <>
450 OCTAVE_API void
451 lu<FloatMatrix>::update (const FloatColumnVector& u,
452  const FloatColumnVector& v)
453 {
454  if (packed ())
455  unpack ();
456 
457  FloatMatrix& l = m_L;
458  FloatMatrix& r = m_a_fact;
459 
460  F77_INT m = to_f77_int (l.rows ());
461  F77_INT n = to_f77_int (r.columns ());
462  F77_INT k = to_f77_int (l.columns ());
463 
464  F77_INT u_nel = to_f77_int (u.numel ());
465  F77_INT v_nel = to_f77_int (v.numel ());
466 
467  if (u_nel != m || v_nel != n)
468  (*current_liboctave_error_handler) ("luupdate: dimensions mismatch");
469 
470  FloatColumnVector utmp = u;
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);
504  F77_XFCN (slu1up, SLU1UP, (m, n, l.fortran_vec (),
505  m, r.fortran_vec (), k,
506  utmp.fortran_vec (), vtmp.fortran_vec ()));
507  }
508 }
509 
510 template <>
511 OCTAVE_API void
512 lu<FloatMatrix>::update_piv (const FloatColumnVector& u,
513  const FloatColumnVector& v)
514 {
515  if (packed ())
516  unpack ();
517 
518  FloatMatrix& l = m_L;
519  FloatMatrix& r = m_a_fact;
520 
521  F77_INT m = to_f77_int (l.rows ());
522  F77_INT n = to_f77_int (r.columns ());
523  F77_INT k = to_f77_int (l.columns ());
524 
525  F77_INT u_nel = to_f77_int (u.numel ());
526  F77_INT v_nel = to_f77_int (v.numel ());
527 
528  if (u_nel != m || v_nel != n)
529  (*current_liboctave_error_handler) ("luupdate: dimensions mismatch");
530 
531  FloatColumnVector utmp = u;
532  FloatColumnVector vtmp = v;
533  OCTAVE_LOCAL_BUFFER (float, w, m);
534  for (F77_INT i = 0; i < m; i++) m_ipvt(i) += 1; // increment
535  F77_XFCN (slup1up, SLUP1UP, (m, n, l.fortran_vec (),
536  m, r.fortran_vec (), k,
537  m_ipvt.fortran_vec (),
538  utmp.data (), vtmp.data (), w));
539  for (F77_INT i = 0; i < m; i++) m_ipvt(i) -= 1; // decrement
540 }
541 
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 
579 #endif
580 
581 template <>
582 OCTAVE_API
583 lu<ComplexMatrix>::lu (const ComplexMatrix& a)
584 {
585  F77_INT a_nr = to_f77_int (a.rows ());
586  F77_INT a_nc = to_f77_int (a.columns ());
587  F77_INT mn = (a_nr < a_nc ? a_nr : a_nc);
588 
589  m_ipvt.resize (dim_vector (mn, 1));
590  F77_INT *pipvt = m_ipvt.fortran_vec ();
591 
592  m_a_fact = a;
593  Complex *tmp_data = m_a_fact.fortran_vec ();
594 
595  F77_INT info = 0;
596 
597  F77_XFCN (zgetrf, ZGETRF, (a_nr, a_nc, F77_DBLE_CMPLX_ARG (tmp_data),
598  a_nr, pipvt, info));
599 
600  for (F77_INT i = 0; i < mn; i++)
601  pipvt[i] -= 1;
602 }
603 
604 #if defined (HAVE_QRUPDATE_LUU)
605 
606 template <>
607 OCTAVE_API void
608 lu<ComplexMatrix>::update (const ComplexColumnVector& u,
609  const ComplexColumnVector& v)
610 {
611  if (packed ())
612  unpack ();
613 
614  ComplexMatrix& l = m_L;
615  ComplexMatrix& r = m_a_fact;
616 
617  F77_INT m = to_f77_int (l.rows ());
618  F77_INT n = to_f77_int (r.columns ());
619  F77_INT k = to_f77_int (l.columns ());
620 
621  F77_INT u_nel = to_f77_int (u.numel ());
622  F77_INT v_nel = to_f77_int (v.numel ());
623 
624  if (u_nel != m || v_nel != n)
625  (*current_liboctave_error_handler) ("luupdate: dimensions mismatch");
626 
627  ComplexColumnVector utmp = u;
628  ComplexColumnVector vtmp = v;
629  F77_XFCN (zlu1up, ZLU1UP, (m, n, F77_DBLE_CMPLX_ARG (l.fortran_vec ()), m,
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,
666  F77_DBLE_CMPLX_ARG (utmp.fortran_vec ()),
667  F77_DBLE_CMPLX_ARG (vtmp.fortran_vec ())));
668  }
669 }
670 
671 template <>
672 OCTAVE_API void
673 lu<ComplexMatrix>::update_piv (const ComplexColumnVector& u,
674  const ComplexColumnVector& v)
675 {
676  if (packed ())
677  unpack ();
678 
679  ComplexMatrix& l = m_L;
680  ComplexMatrix& r = m_a_fact;
681 
682  F77_INT m = to_f77_int (l.rows ());
683  F77_INT n = to_f77_int (r.columns ());
684  F77_INT k = to_f77_int (l.columns ());
685 
686  F77_INT u_nel = to_f77_int (u.numel ());
687  F77_INT v_nel = to_f77_int (v.numel ());
688 
689  if (u_nel != m || v_nel != n)
690  (*current_liboctave_error_handler) ("luupdate: dimensions mismatch");
691 
692  ComplexColumnVector utmp = u;
693  ComplexColumnVector vtmp = v;
694  OCTAVE_LOCAL_BUFFER (Complex, w, m);
695  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 ()),
697  m, F77_DBLE_CMPLX_ARG (r.fortran_vec ()), k,
698  m_ipvt.fortran_vec (),
699  F77_CONST_DBLE_CMPLX_ARG (utmp.data ()),
700  F77_CONST_DBLE_CMPLX_ARG (vtmp.data ()),
701  F77_DBLE_CMPLX_ARG (w)));
702  for (F77_INT i = 0; i < m; i++) m_ipvt(i) -= 1; // decrement
703 }
704 
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 
747 #endif
748 
749 template <>
750 OCTAVE_API
751 lu<FloatComplexMatrix>::lu (const FloatComplexMatrix& a)
752 {
753  F77_INT a_nr = to_f77_int (a.rows ());
754  F77_INT a_nc = to_f77_int (a.columns ());
755  F77_INT mn = (a_nr < a_nc ? a_nr : a_nc);
756 
757  m_ipvt.resize (dim_vector (mn, 1));
758  F77_INT *pipvt = m_ipvt.fortran_vec ();
759 
760  m_a_fact = a;
761  FloatComplex *tmp_data = m_a_fact.fortran_vec ();
762 
763  F77_INT info = 0;
764 
765  F77_XFCN (cgetrf, CGETRF, (a_nr, a_nc, F77_CMPLX_ARG (tmp_data), a_nr,
766  pipvt, info));
767 
768  for (F77_INT i = 0; i < mn; i++)
769  pipvt[i] -= 1;
770 }
771 
772 #if defined (HAVE_QRUPDATE_LUU)
773 
774 template <>
775 OCTAVE_API void
776 lu<FloatComplexMatrix>::update (const FloatComplexColumnVector& u,
777  const FloatComplexColumnVector& v)
778 {
779  if (packed ())
780  unpack ();
781 
782  FloatComplexMatrix& l = m_L;
783  FloatComplexMatrix& r = m_a_fact;
784 
785  F77_INT m = to_f77_int (l.rows ());
786  F77_INT n = to_f77_int (r.columns ());
787  F77_INT k = to_f77_int (l.columns ());
788 
789  F77_INT u_nel = to_f77_int (u.numel ());
790  F77_INT v_nel = to_f77_int (v.numel ());
791 
792  if (u_nel != m || v_nel != n)
793  (*current_liboctave_error_handler) ("luupdate: dimensions mismatch");
794 
795  FloatComplexColumnVector utmp = u;
796  FloatComplexColumnVector vtmp = v;
797  F77_XFCN (clu1up, CLU1UP, (m, n, F77_CMPLX_ARG (l.fortran_vec ()), m,
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,
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 ()),
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,
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 
912 #endif
913 
914 // Instantiations we need.
915 
916 template class lu<Matrix>;
917 
918 template class lu<FloatMatrix>;
919 
920 template class lu<ComplexMatrix>;
921 
922 template class lu<FloatComplexMatrix>;
923 
924 OCTAVE_END_NAMESPACE(math)
925 OCTAVE_END_NAMESPACE(octave)
min
charNDArray min(char d, const charNDArray &m)
Definition: chNDArray.cc:207
Array::fortran_vec
T * fortran_vec()
Size of the specified dimension.
Definition: Array-base.cc:1764
Array::rows
octave_idx_type rows() const
Definition: Array.h:459
Array::data
const T * data() const
Size of the specified dimension.
Definition: Array.h:663
Array::columns
octave_idx_type columns() const
Definition: Array.h:471
Array::xelem
T & xelem(octave_idx_type n)
Size of the specified dimension.
Definition: Array.h:524
Array::numel
octave_idx_type numel() const
Number of elements in the array.
Definition: Array.h:414
ComplexMatrix::column
ComplexColumnVector column(octave_idx_type i) const
Definition: CMatrix.cc:708
FloatComplexMatrix::column
FloatComplexColumnVector column(octave_idx_type i) const
Definition: fCMatrix.cc:711
FloatMatrix::column
FloatColumnVector column(octave_idx_type i) const
Definition: fMatrix.cc:428
Matrix
Definition: dMatrix.h:42
Matrix::column
ColumnVector column(octave_idx_type i) const
Definition: dMatrix.cc:422
dim_vector
Vector representing the dimensions (size) of an Array.
Definition: dim-vector.h:94
lu
Definition: lu.h:42
lu::packed
bool packed() const
Definition: lu.cc:68
lu::lu
lu()
Definition: lu.h:48
lu::update_piv
void update_piv(const VT &u, const VT &v)
lu::U
T U() const
Definition: lu.cc:123
lu::unpack
void unpack()
Definition: lu.cc:75
lu::P_vec
ColumnVector P_vec() const
Definition: lu.cc:196
lu::regular
bool regular() const
Definition: lu.cc:212
lu::update
void update(const VT &u, const VT &v)
lu::ELT_T
T::element_type ELT_T
Definition: lu.h:46
lu::P
PermMatrix P() const
Definition: lu.cc:189
lu::getp
Array< octave_idx_type > getp() const
Definition: lu.cc:158
lu::L
T L() const
Definition: lu.cc:96
lu::Y
T Y() const
Definition: lu.cc:147
octave_idx_type
OCTAVE_BEGIN_NAMESPACE
OCTAVE_BEGIN_NAMESPACE(octave) static octave_value daspk_fcn
F77_CONST_CMPLX_ARG
#define F77_CONST_CMPLX_ARG(x)
Definition: f77-fcn.h:313
F77_DBLE_CMPLX_ARG
#define F77_DBLE_CMPLX_ARG(x)
Definition: f77-fcn.h:316
F77_CMPLX_ARG
#define F77_CMPLX_ARG(x)
Definition: f77-fcn.h:310
F77_XFCN
#define F77_XFCN(f, F, args)
Definition: f77-fcn.h:45
F77_INT
octave_f77_int_type F77_INT
Definition: f77-fcn.h:306
F77_CONST_DBLE_CMPLX_ARG
#define F77_CONST_DBLE_CMPLX_ARG(x)
Definition: f77-fcn.h:319
OCTAVE_API
#define OCTAVE_API
Definition: main.cc:55
m
T octave_idx_type m
Definition: mx-inlines.cc:781
n
octave_idx_type n
Definition: mx-inlines.cc:761
r
T * r
Definition: mx-inlines.cc:781
Faddeeva::w
std::complex< double > w(std::complex< double > z, double relerr=0)
Complex
std::complex< double > Complex
Definition: oct-cmplx.h:33
FloatComplex
std::complex< float > FloatComplex
Definition: oct-cmplx.h:34
OCTAVE_LOCAL_BUFFER
#define OCTAVE_LOCAL_BUFFER(T, buf, size)
Definition: oct-locbuf.h:44