GNU Octave  9.1.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
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