GNU Octave  6.2.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-2021 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 namespace octave
48 {
49  namespace math
50  {
51  // FIXME: PermMatrix::col_perm_vec returns Array<octave_idx_type>
52  // but 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  : a_fact (u), l_fact (l), 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 (void) const
69  {
70  return l_fact.dims () == dim_vector ();
71  }
72 
73  template <typename T>
74  void
75  lu<T>::unpack (void)
76  {
77  if (packed ())
78  {
79  l_fact = L ();
80  a_fact = U (); // FIXME: sub-optimal
81 
82  // FIXME: getp returns Array<octave_idx_type> but ipvt is
83  // Array<octave_f77_int_type>. However, getp produces its
84  // result from a valid 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  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 = a_fact.rows ();
101  octave_idx_type a_nc = 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) = a_fact.xelem (i, j);
113  }
114 
115  return l;
116  }
117  else
118  return l_fact;
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 = a_fact.rows ();
128  octave_idx_type a_nc = 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) = a_fact.xelem (i, j);
137  }
138 
139  return u;
140  }
141  else
142  return 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 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 = 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 < ipvt.numel (); i++)
170  {
171  octave_idx_type k = 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 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 = 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 (a_fact.rows (), a_fact.columns ());
217 
218  for (octave_idx_type i = 0; i < k; i++)
219  {
220  if (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  lu<Matrix>::lu (const Matrix& a)
274  {
275  F77_INT a_nr = to_f77_int (a.rows ());
276  F77_INT a_nc = to_f77_int (a.columns ());
277  F77_INT mn = (a_nr < a_nc ? a_nr : a_nc);
278 
279  ipvt.resize (dim_vector (mn, 1));
280  F77_INT *pipvt = ipvt.fortran_vec ();
281 
282  a_fact = a;
283  double *tmp_data = a_fact.fortran_vec ();
284 
285  F77_INT info = 0;
286 
287  F77_XFCN (dgetrf, DGETRF, (a_nr, a_nc, tmp_data, a_nr, pipvt, info));
288 
289  for (F77_INT i = 0; i < mn; i++)
290  pipvt[i] -= 1;
291  }
292 
293 #if defined (HAVE_QRUPDATE_LUU)
294 
295  template <>
296  void
297  lu<Matrix>::update (const ColumnVector& u, const ColumnVector& v)
298  {
299  if (packed ())
300  unpack ();
301 
302  Matrix& l = l_fact;
303  Matrix& r = a_fact;
304 
305  F77_INT m = to_f77_int (l.rows ());
306  F77_INT n = to_f77_int (r.columns ());
307  F77_INT k = to_f77_int (l.columns ());
308 
309  F77_INT u_nel = to_f77_int (u.numel ());
310  F77_INT v_nel = to_f77_int (v.numel ());
311 
312  if (u_nel != m || v_nel != n)
313  (*current_liboctave_error_handler) ("luupdate: dimensions mismatch");
314 
315  ColumnVector utmp = u;
316  ColumnVector vtmp = v;
317  F77_XFCN (dlu1up, DLU1UP, (m, n, l.fortran_vec (), m, r.fortran_vec (),
318  k, utmp.fortran_vec (), vtmp.fortran_vec ()));
319  }
320 
321  template <>
322  void
323  lu<Matrix>::update (const Matrix& u, const Matrix& v)
324  {
325  if (packed ())
326  unpack ();
327 
328  Matrix& l = l_fact;
329  Matrix& r = a_fact;
330 
331  F77_INT m = to_f77_int (l.rows ());
332  F77_INT n = to_f77_int (r.columns ());
333  F77_INT k = to_f77_int (l.columns ());
334 
335  F77_INT u_nr = to_f77_int (u.rows ());
336  F77_INT u_nc = to_f77_int (u.columns ());
337 
338  F77_INT v_nr = to_f77_int (v.rows ());
339  F77_INT v_nc = to_f77_int (v.columns ());
340 
341  if (u_nr != m || v_nr != n || u_nc != v_nc)
342  (*current_liboctave_error_handler) ("luupdate: dimensions mismatch");
343 
344  for (volatile F77_INT i = 0; i < u_nc; i++)
345  {
346  ColumnVector utmp = u.column (i);
347  ColumnVector vtmp = v.column (i);
348  F77_XFCN (dlu1up, DLU1UP, (m, n, l.fortran_vec (),
349  m, r.fortran_vec (), k,
350  utmp.fortran_vec (), vtmp.fortran_vec ()));
351  }
352  }
353 
354  template <>
355  void
356  lu<Matrix>::update_piv (const ColumnVector& u, const ColumnVector& v)
357  {
358  if (packed ())
359  unpack ();
360 
361  Matrix& l = l_fact;
362  Matrix& r = a_fact;
363 
364  F77_INT m = to_f77_int (l.rows ());
365  F77_INT n = to_f77_int (r.columns ());
366  F77_INT k = to_f77_int (l.columns ());
367 
368  F77_INT u_nel = to_f77_int (u.numel ());
369  F77_INT v_nel = to_f77_int (v.numel ());
370 
371  if (u_nel != m || v_nel != n)
372  (*current_liboctave_error_handler) ("luupdate: dimensions mismatch");
373 
374  ColumnVector utmp = u;
375  ColumnVector vtmp = v;
376  OCTAVE_LOCAL_BUFFER (double, w, m);
377  for (F77_INT i = 0; i < m; i++) ipvt(i) += 1; // increment
378  F77_XFCN (dlup1up, DLUP1UP, (m, n, l.fortran_vec (),
379  m, r.fortran_vec (), k,
380  ipvt.fortran_vec (),
381  utmp.data (), vtmp.data (), w));
382  for (F77_INT i = 0; i < m; i++) ipvt(i) -= 1; // decrement
383  }
384 
385  template <>
386  void
387  lu<Matrix>::update_piv (const Matrix& u, const Matrix& v)
388  {
389  if (packed ())
390  unpack ();
391 
392  Matrix& l = l_fact;
393  Matrix& r = a_fact;
394 
395  F77_INT m = to_f77_int (l.rows ());
396  F77_INT n = to_f77_int (r.columns ());
397  F77_INT k = to_f77_int (l.columns ());
398 
399  F77_INT u_nr = to_f77_int (u.rows ());
400  F77_INT u_nc = to_f77_int (u.columns ());
401 
402  F77_INT v_nr = to_f77_int (v.rows ());
403  F77_INT v_nc = to_f77_int (v.columns ());
404 
405  if (u_nr != m || v_nr != n || u_nc != v_nc)
406  (*current_liboctave_error_handler) ("luupdate: dimensions mismatch");
407 
408  OCTAVE_LOCAL_BUFFER (double, w, m);
409  for (F77_INT i = 0; i < m; i++) ipvt(i) += 1; // increment
410  for (volatile F77_INT i = 0; i < u_nc; i++)
411  {
412  ColumnVector utmp = u.column (i);
413  ColumnVector vtmp = v.column (i);
414  F77_XFCN (dlup1up, DLUP1UP, (m, n, l.fortran_vec (),
415  m, r.fortran_vec (), k,
416  ipvt.fortran_vec (),
417  utmp.data (), vtmp.data (), w));
418  }
419  for (F77_INT i = 0; i < m; i++) ipvt(i) -= 1; // decrement
420  }
421 
422 #endif
423 
424  template <>
425  lu<FloatMatrix>::lu (const FloatMatrix& a)
426  {
427  F77_INT a_nr = to_f77_int (a.rows ());
428  F77_INT a_nc = to_f77_int (a.columns ());
429  F77_INT mn = (a_nr < a_nc ? a_nr : a_nc);
430 
431  ipvt.resize (dim_vector (mn, 1));
432  F77_INT *pipvt = ipvt.fortran_vec ();
433 
434  a_fact = a;
435  float *tmp_data = a_fact.fortran_vec ();
436 
437  F77_INT info = 0;
438 
439  F77_XFCN (sgetrf, SGETRF, (a_nr, a_nc, tmp_data, a_nr, pipvt, info));
440 
441  for (F77_INT i = 0; i < mn; i++)
442  pipvt[i] -= 1;
443  }
444 
445 #if defined (HAVE_QRUPDATE_LUU)
446 
447  template <>
448  void
449  lu<FloatMatrix>::update (const FloatColumnVector& u,
450  const FloatColumnVector& v)
451  {
452  if (packed ())
453  unpack ();
454 
455  FloatMatrix& l = l_fact;
456  FloatMatrix& r = a_fact;
457 
458  F77_INT m = to_f77_int (l.rows ());
459  F77_INT n = to_f77_int (r.columns ());
460  F77_INT k = to_f77_int (l.columns ());
461 
462  F77_INT u_nel = to_f77_int (u.numel ());
463  F77_INT v_nel = to_f77_int (v.numel ());
464 
465  if (u_nel != m || v_nel != n)
466  (*current_liboctave_error_handler) ("luupdate: dimensions mismatch");
467 
468  FloatColumnVector utmp = u;
469  FloatColumnVector vtmp = v;
470  F77_XFCN (slu1up, SLU1UP, (m, n, l.fortran_vec (),
471  m, r.fortran_vec (), k,
472  utmp.fortran_vec (), vtmp.fortran_vec ()));
473  }
474 
475  template <>
476  void
477  lu<FloatMatrix>::update (const FloatMatrix& u, const FloatMatrix& v)
478  {
479  if (packed ())
480  unpack ();
481 
482  FloatMatrix& l = l_fact;
483  FloatMatrix& r = a_fact;
484 
485  F77_INT m = to_f77_int (l.rows ());
486  F77_INT n = to_f77_int (r.columns ());
487  F77_INT k = to_f77_int (l.columns ());
488 
489  F77_INT u_nr = to_f77_int (u.rows ());
490  F77_INT u_nc = to_f77_int (u.columns ());
491 
492  F77_INT v_nr = to_f77_int (v.rows ());
493  F77_INT v_nc = to_f77_int (v.columns ());
494 
495  if (u_nr != m || v_nr != n || u_nc != v_nc)
496  (*current_liboctave_error_handler) ("luupdate: dimensions mismatch");
497 
498  for (volatile F77_INT i = 0; i < u_nc; i++)
499  {
500  FloatColumnVector utmp = u.column (i);
501  FloatColumnVector vtmp = v.column (i);
502  F77_XFCN (slu1up, SLU1UP, (m, n, l.fortran_vec (),
503  m, r.fortran_vec (), k,
504  utmp.fortran_vec (), vtmp.fortran_vec ()));
505  }
506  }
507 
508  template <>
509  void
510  lu<FloatMatrix>::update_piv (const FloatColumnVector& u,
511  const FloatColumnVector& v)
512  {
513  if (packed ())
514  unpack ();
515 
516  FloatMatrix& l = l_fact;
517  FloatMatrix& r = a_fact;
518 
519  F77_INT m = to_f77_int (l.rows ());
520  F77_INT n = to_f77_int (r.columns ());
521  F77_INT k = to_f77_int (l.columns ());
522 
523  F77_INT u_nel = to_f77_int (u.numel ());
524  F77_INT v_nel = to_f77_int (v.numel ());
525 
526  if (u_nel != m || v_nel != n)
527  (*current_liboctave_error_handler) ("luupdate: dimensions mismatch");
528 
529  FloatColumnVector utmp = u;
530  FloatColumnVector vtmp = v;
531  OCTAVE_LOCAL_BUFFER (float, w, m);
532  for (F77_INT i = 0; i < m; i++) ipvt(i) += 1; // increment
533  F77_XFCN (slup1up, SLUP1UP, (m, n, l.fortran_vec (),
534  m, r.fortran_vec (), k,
535  ipvt.fortran_vec (),
536  utmp.data (), vtmp.data (), w));
537  for (F77_INT i = 0; i < m; i++) ipvt(i) -= 1; // decrement
538  }
539 
540  template <>
541  void
542  lu<FloatMatrix>::update_piv (const FloatMatrix& u, const FloatMatrix& v)
543  {
544  if (packed ())
545  unpack ();
546 
547  FloatMatrix& l = l_fact;
548  FloatMatrix& r = a_fact;
549 
550  F77_INT m = to_f77_int (l.rows ());
551  F77_INT n = to_f77_int (r.columns ());
552  F77_INT k = to_f77_int (l.columns ());
553 
554  F77_INT u_nr = to_f77_int (u.rows ());
555  F77_INT u_nc = to_f77_int (u.columns ());
556 
557  F77_INT v_nr = to_f77_int (v.rows ());
558  F77_INT v_nc = to_f77_int (v.columns ());
559 
560  if (u_nr != m || v_nr != n || u_nc != v_nc)
561  (*current_liboctave_error_handler) ("luupdate: dimensions mismatch");
562 
563  OCTAVE_LOCAL_BUFFER (float, w, m);
564  for (F77_INT i = 0; i < m; i++) ipvt(i) += 1; // increment
565  for (volatile F77_INT i = 0; i < u_nc; i++)
566  {
567  FloatColumnVector utmp = u.column (i);
568  FloatColumnVector vtmp = v.column (i);
569  F77_XFCN (slup1up, SLUP1UP, (m, n, l.fortran_vec (),
570  m, r.fortran_vec (), k,
571  ipvt.fortran_vec (),
572  utmp.data (), vtmp.data (), w));
573  }
574  for (F77_INT i = 0; i < m; i++) ipvt(i) -= 1; // decrement
575  }
576 
577 #endif
578 
579  template <>
580  lu<ComplexMatrix>::lu (const ComplexMatrix& a)
581  {
582  F77_INT a_nr = to_f77_int (a.rows ());
583  F77_INT a_nc = to_f77_int (a.columns ());
584  F77_INT mn = (a_nr < a_nc ? a_nr : a_nc);
585 
586  ipvt.resize (dim_vector (mn, 1));
587  F77_INT *pipvt = ipvt.fortran_vec ();
588 
589  a_fact = a;
590  Complex *tmp_data = a_fact.fortran_vec ();
591 
592  F77_INT info = 0;
593 
594  F77_XFCN (zgetrf, ZGETRF, (a_nr, a_nc, F77_DBLE_CMPLX_ARG (tmp_data),
595  a_nr, pipvt, info));
596 
597  for (F77_INT i = 0; i < mn; i++)
598  pipvt[i] -= 1;
599  }
600 
601 #if defined (HAVE_QRUPDATE_LUU)
602 
603  template <>
604  void
605  lu<ComplexMatrix>::update (const ComplexColumnVector& u,
606  const ComplexColumnVector& v)
607  {
608  if (packed ())
609  unpack ();
610 
611  ComplexMatrix& l = l_fact;
612  ComplexMatrix& r = a_fact;
613 
614  F77_INT m = to_f77_int (l.rows ());
615  F77_INT n = to_f77_int (r.columns ());
616  F77_INT k = to_f77_int (l.columns ());
617 
618  F77_INT u_nel = to_f77_int (u.numel ());
619  F77_INT v_nel = to_f77_int (v.numel ());
620 
621  if (u_nel != m || v_nel != n)
622  (*current_liboctave_error_handler) ("luupdate: dimensions mismatch");
623 
624  ComplexColumnVector utmp = u;
625  ComplexColumnVector vtmp = v;
626  F77_XFCN (zlu1up, ZLU1UP, (m, n, F77_DBLE_CMPLX_ARG (l.fortran_vec ()), m,
627  F77_DBLE_CMPLX_ARG (r.fortran_vec ()), k,
628  F77_DBLE_CMPLX_ARG (utmp.fortran_vec ()),
629  F77_DBLE_CMPLX_ARG (vtmp.fortran_vec ())));
630  }
631 
632  template <>
633  void
634  lu<ComplexMatrix>::update (const ComplexMatrix& u, const ComplexMatrix& v)
635  {
636  if (packed ())
637  unpack ();
638 
639  ComplexMatrix& l = l_fact;
640  ComplexMatrix& r = a_fact;
641 
642  F77_INT m = to_f77_int (l.rows ());
643  F77_INT n = to_f77_int (r.columns ());
644  F77_INT k = to_f77_int (l.columns ());
645 
646  F77_INT u_nr = to_f77_int (u.rows ());
647  F77_INT u_nc = to_f77_int (u.columns ());
648 
649  F77_INT v_nr = to_f77_int (v.rows ());
650  F77_INT v_nc = to_f77_int (v.columns ());
651 
652  if (u_nr != m || v_nr != n || u_nc != v_nc)
653  (*current_liboctave_error_handler) ("luupdate: dimensions mismatch");
654 
655  for (volatile F77_INT i = 0; i < u_nc; i++)
656  {
657  ComplexColumnVector utmp = u.column (i);
658  ComplexColumnVector vtmp = v.column (i);
659  F77_XFCN (zlu1up, ZLU1UP, (m, n,
660  F77_DBLE_CMPLX_ARG (l.fortran_vec ()),
661  m, F77_DBLE_CMPLX_ARG (r.fortran_vec ()),
662  k,
663  F77_DBLE_CMPLX_ARG (utmp.fortran_vec ()),
664  F77_DBLE_CMPLX_ARG (vtmp.fortran_vec ())));
665  }
666  }
667 
668  template <>
669  void
670  lu<ComplexMatrix>::update_piv (const ComplexColumnVector& u,
671  const ComplexColumnVector& v)
672  {
673  if (packed ())
674  unpack ();
675 
676  ComplexMatrix& l = l_fact;
677  ComplexMatrix& r = a_fact;
678 
679  F77_INT m = to_f77_int (l.rows ());
680  F77_INT n = to_f77_int (r.columns ());
681  F77_INT k = to_f77_int (l.columns ());
682 
683  F77_INT u_nel = to_f77_int (u.numel ());
684  F77_INT v_nel = to_f77_int (v.numel ());
685 
686  if (u_nel != m || v_nel != n)
687  (*current_liboctave_error_handler) ("luupdate: dimensions mismatch");
688 
689  ComplexColumnVector utmp = u;
690  ComplexColumnVector vtmp = v;
691  OCTAVE_LOCAL_BUFFER (Complex, w, m);
692  for (F77_INT i = 0; i < m; i++) ipvt(i) += 1; // increment
693  F77_XFCN (zlup1up, ZLUP1UP, (m, n, F77_DBLE_CMPLX_ARG (l.fortran_vec ()),
694  m, F77_DBLE_CMPLX_ARG (r.fortran_vec ()), k,
695  ipvt.fortran_vec (),
696  F77_CONST_DBLE_CMPLX_ARG (utmp.data ()),
697  F77_CONST_DBLE_CMPLX_ARG (vtmp.data ()),
698  F77_DBLE_CMPLX_ARG (w)));
699  for (F77_INT i = 0; i < m; i++) ipvt(i) -= 1; // decrement
700  }
701 
702  template <>
703  void
704  lu<ComplexMatrix>::update_piv (const ComplexMatrix& u,
705  const ComplexMatrix& v)
706  {
707  if (packed ())
708  unpack ();
709 
710  ComplexMatrix& l = l_fact;
711  ComplexMatrix& r = a_fact;
712 
713  F77_INT m = to_f77_int (l.rows ());
714  F77_INT n = to_f77_int (r.columns ());
715  F77_INT k = to_f77_int (l.columns ());
716 
717  F77_INT u_nr = to_f77_int (u.rows ());
718  F77_INT u_nc = to_f77_int (u.columns ());
719 
720  F77_INT v_nr = to_f77_int (v.rows ());
721  F77_INT v_nc = to_f77_int (v.columns ());
722 
723  if (u_nr != m || v_nr != n || u_nc != v_nc)
724  (*current_liboctave_error_handler) ("luupdate: dimensions mismatch");
725 
726  OCTAVE_LOCAL_BUFFER (Complex, w, m);
727  for (F77_INT i = 0; i < m; i++) ipvt(i) += 1; // increment
728  for (volatile F77_INT i = 0; i < u_nc; i++)
729  {
730  ComplexColumnVector utmp = u.column (i);
731  ComplexColumnVector vtmp = v.column (i);
732  F77_XFCN (zlup1up, ZLUP1UP, (m, n,
733  F77_DBLE_CMPLX_ARG (l.fortran_vec ()),
734  m,
735  F77_DBLE_CMPLX_ARG (r.fortran_vec ()),
736  k, ipvt.fortran_vec (),
737  F77_CONST_DBLE_CMPLX_ARG (utmp.data ()),
738  F77_CONST_DBLE_CMPLX_ARG (vtmp.data ()),
739  F77_DBLE_CMPLX_ARG (w)));
740  }
741  for (F77_INT i = 0; i < m; i++) ipvt(i) -= 1; // decrement
742  }
743 
744 #endif
745 
746  template <>
747  lu<FloatComplexMatrix>::lu (const FloatComplexMatrix& a)
748  {
749  F77_INT a_nr = to_f77_int (a.rows ());
750  F77_INT a_nc = to_f77_int (a.columns ());
751  F77_INT mn = (a_nr < a_nc ? a_nr : a_nc);
752 
753  ipvt.resize (dim_vector (mn, 1));
754  F77_INT *pipvt = ipvt.fortran_vec ();
755 
756  a_fact = a;
757  FloatComplex *tmp_data = a_fact.fortran_vec ();
758 
759  F77_INT info = 0;
760 
761  F77_XFCN (cgetrf, CGETRF, (a_nr, a_nc, F77_CMPLX_ARG (tmp_data), a_nr,
762  pipvt, info));
763 
764  for (F77_INT i = 0; i < mn; i++)
765  pipvt[i] -= 1;
766  }
767 
768 #if defined (HAVE_QRUPDATE_LUU)
769 
770  template <>
771  void
772  lu<FloatComplexMatrix>::update (const FloatComplexColumnVector& u,
773  const FloatComplexColumnVector& v)
774  {
775  if (packed ())
776  unpack ();
777 
778  FloatComplexMatrix& l = l_fact;
779  FloatComplexMatrix& r = a_fact;
780 
781  F77_INT m = to_f77_int (l.rows ());
782  F77_INT n = to_f77_int (r.columns ());
783  F77_INT k = to_f77_int (l.columns ());
784 
785  F77_INT u_nel = to_f77_int (u.numel ());
786  F77_INT v_nel = to_f77_int (v.numel ());
787 
788  if (u_nel != m || v_nel != n)
789  (*current_liboctave_error_handler) ("luupdate: dimensions mismatch");
790 
791  FloatComplexColumnVector utmp = u;
792  FloatComplexColumnVector vtmp = v;
793  F77_XFCN (clu1up, CLU1UP, (m, n, F77_CMPLX_ARG (l.fortran_vec ()), m,
794  F77_CMPLX_ARG (r.fortran_vec ()), k,
795  F77_CMPLX_ARG (utmp.fortran_vec ()),
796  F77_CMPLX_ARG (vtmp.fortran_vec ())));
797  }
798 
799  template <>
800  void
801  lu<FloatComplexMatrix>::update (const FloatComplexMatrix& u,
802  const FloatComplexMatrix& v)
803  {
804  if (packed ())
805  unpack ();
806 
807  FloatComplexMatrix& l = l_fact;
808  FloatComplexMatrix& r = a_fact;
809 
810  F77_INT m = to_f77_int (l.rows ());
811  F77_INT n = to_f77_int (r.columns ());
812  F77_INT k = to_f77_int (l.columns ());
813 
814  F77_INT u_nr = to_f77_int (u.rows ());
815  F77_INT u_nc = to_f77_int (u.columns ());
816 
817  F77_INT v_nr = to_f77_int (v.rows ());
818  F77_INT v_nc = to_f77_int (v.columns ());
819 
820  if (u_nr != m || v_nr != n || u_nc != v_nc)
821  (*current_liboctave_error_handler) ("luupdate: dimensions mismatch");
822 
823  for (volatile F77_INT i = 0; i < u_nc; i++)
824  {
825  FloatComplexColumnVector utmp = u.column (i);
826  FloatComplexColumnVector vtmp = v.column (i);
827  F77_XFCN (clu1up, CLU1UP, (m, n, F77_CMPLX_ARG (l.fortran_vec ()),
828  m, F77_CMPLX_ARG (r.fortran_vec ()), k,
829  F77_CMPLX_ARG (utmp.fortran_vec ()),
830  F77_CMPLX_ARG (vtmp.fortran_vec ())));
831  }
832  }
833 
834  template <>
835  void
836  lu<FloatComplexMatrix>::update_piv (const FloatComplexColumnVector& u,
837  const FloatComplexColumnVector& v)
838  {
839  if (packed ())
840  unpack ();
841 
842  FloatComplexMatrix& l = l_fact;
843  FloatComplexMatrix& r = a_fact;
844 
845  F77_INT m = to_f77_int (l.rows ());
846  F77_INT n = to_f77_int (r.columns ());
847  F77_INT k = to_f77_int (l.columns ());
848 
849  F77_INT u_nel = to_f77_int (u.numel ());
850  F77_INT v_nel = to_f77_int (v.numel ());
851 
852  if (u_nel != m || v_nel != n)
853  (*current_liboctave_error_handler) ("luupdate: dimensions mismatch");
854 
855  FloatComplexColumnVector utmp = u;
856  FloatComplexColumnVector vtmp = v;
857  OCTAVE_LOCAL_BUFFER (FloatComplex, w, m);
858  for (F77_INT i = 0; i < m; i++) ipvt(i) += 1; // increment
859  F77_XFCN (clup1up, CLUP1UP, (m, n, F77_CMPLX_ARG (l.fortran_vec ()),
860  m, F77_CMPLX_ARG (r.fortran_vec ()), k,
861  ipvt.fortran_vec (),
862  F77_CONST_CMPLX_ARG (utmp.data ()),
863  F77_CONST_CMPLX_ARG (vtmp.data ()),
864  F77_CMPLX_ARG (w)));
865  for (F77_INT i = 0; i < m; i++) ipvt(i) -= 1; // decrement
866  }
867 
868  template <>
869  void
870  lu<FloatComplexMatrix>::update_piv (const FloatComplexMatrix& u,
871  const FloatComplexMatrix& v)
872  {
873  if (packed ())
874  unpack ();
875 
876  FloatComplexMatrix& l = l_fact;
877  FloatComplexMatrix& r = a_fact;
878 
879  F77_INT m = to_f77_int (l.rows ());
880  F77_INT n = to_f77_int (r.columns ());
881  F77_INT k = to_f77_int (l.columns ());
882 
883  F77_INT u_nr = to_f77_int (u.rows ());
884  F77_INT u_nc = to_f77_int (u.columns ());
885 
886  F77_INT v_nr = to_f77_int (v.rows ());
887  F77_INT v_nc = to_f77_int (v.columns ());
888 
889  if (u_nr != m || v_nr != n || u_nc != v_nc)
890  (*current_liboctave_error_handler) ("luupdate: dimensions mismatch");
891 
892  OCTAVE_LOCAL_BUFFER (FloatComplex, w, m);
893  for (F77_INT i = 0; i < m; i++) ipvt(i) += 1; // increment
894  for (volatile F77_INT i = 0; i < u_nc; i++)
895  {
896  FloatComplexColumnVector utmp = u.column (i);
897  FloatComplexColumnVector vtmp = v.column (i);
898  F77_XFCN (clup1up, CLUP1UP, (m, n, F77_CMPLX_ARG (l.fortran_vec ()),
899  m, F77_CMPLX_ARG (r.fortran_vec ()), k,
900  ipvt.fortran_vec (),
901  F77_CONST_CMPLX_ARG (utmp.data ()),
902  F77_CONST_CMPLX_ARG (vtmp.data ()),
903  F77_CMPLX_ARG (w)));
904  }
905  for (F77_INT i = 0; i < m; i++) ipvt(i) -= 1; // decrement
906  }
907 
908 #endif
909 
910  // Instantiations we need.
911 
912  template class lu<Matrix>;
913 
914  template class lu<FloatMatrix>;
915 
916  template class lu<ComplexMatrix>;
917 
918  template class lu<FloatComplexMatrix>;
919  }
920 }
min
charNDArray min(char d, const charNDArray &m)
Definition: chNDArray.cc:207
Array::columns
octave_idx_type columns(void) const
Definition: Array.h:424
Array::xelem
T & xelem(octave_idx_type n)
Size of the specified dimension.
Definition: Array.h:469
Array::numel
octave_idx_type numel(void) const
Number of elements in the array.
Definition: Array.h:377
Array::data
const T * data(void) const
Size of the specified dimension.
Definition: Array.h:581
Array::rows
octave_idx_type rows(void) const
Definition: Array.h:415
Array::fortran_vec
const T * fortran_vec(void) const
Size of the specified dimension.
Definition: Array.h:583
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:95
lu
Definition: mx-defs.h:80
octave::math::lu::lu
lu(void)
Definition: lu.h:49
octave::math::lu::U
T U(void) const
Definition: lu.cc:123
octave::math::lu::regular
bool regular(void) const
Definition: lu.cc:212
octave::math::lu::P
PermMatrix P(void) const
Definition: lu.cc:189
octave::math::lu::update
void update(const VT &u, const VT &v)
octave::math::lu::Y
T Y(void) const
Definition: lu.cc:147
octave::math::lu::unpack
void unpack(void)
Definition: lu.cc:75
octave::math::lu::packed
bool packed(void) const
Definition: lu.cc:68
octave::math::lu::ELT_T
T::element_type ELT_T
Definition: lu.h:47
octave::math::lu::update_piv
void update_piv(const VT &u, const VT &v)
octave::math::lu::P_vec
ColumnVector P_vec(void) const
Definition: lu.cc:196
octave::math::lu::getp
Array< octave_idx_type > getp(void) const
Definition: lu.cc:158
octave::math::lu::L
T L(void) const
Definition: lu.cc:96
octave_idx_type
F77_CONST_CMPLX_ARG
#define F77_CONST_CMPLX_ARG(x)
Definition: f77-fcn.h:312
F77_DBLE_CMPLX_ARG
#define F77_DBLE_CMPLX_ARG(x)
Definition: f77-fcn.h:315
F77_CMPLX_ARG
#define F77_CMPLX_ARG(x)
Definition: f77-fcn.h:309
F77_XFCN
#define F77_XFCN(f, F, args)
Definition: f77-fcn.h:44
F77_INT
octave_f77_int_type F77_INT
Definition: f77-fcn.h:305
F77_CONST_DBLE_CMPLX_ARG
#define F77_CONST_DBLE_CMPLX_ARG(x)
Definition: f77-fcn.h:318
m
T octave_idx_type m
Definition: mx-inlines.cc:773
n
octave_idx_type n
Definition: mx-inlines.cc:753
r
T * r
Definition: mx-inlines.cc:773
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
retval
octave_value::octave_value(const Array< char > &chm, char type) return retval
Definition: ov.cc:811
transpose
static void transpose(octave_idx_type N, const octave_idx_type *ridx, const octave_idx_type *cidx, octave_idx_type *ridx2, octave_idx_type *cidx2)
Definition: symrcm.cc:389