GNU Octave  4.2.1
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
lu.cc
Go to the documentation of this file.
1 /*
2 
3 Copyright (C) 1996-2017 John W. Eaton
4 Copyright (C) 2009 VZLU Prague
5 
6 This file is part of Octave.
7 
8 Octave is free software; you can redistribute it and/or modify it
9 under the terms of the GNU General Public License as published by the
10 Free Software Foundation; either version 3 of the License, or (at your
11 option) any later version.
12 
13 Octave is distributed in the hope that it will be useful, but WITHOUT
14 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
15 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
16 for more details.
17 
18 You should have received a copy of the GNU General Public License
19 along with Octave; see the file COPYING. If not, see
20 <http://www.gnu.org/licenses/>.
21 
22 */
23 
24 #if defined (HAVE_CONFIG_H)
25 # include "config.h"
26 #endif
27 
28 #include "CColVector.h"
29 #include "CMatrix.h"
30 #include "dColVector.h"
31 #include "dMatrix.h"
32 #include "fCColVector.h"
33 #include "fCMatrix.h"
34 #include "fColVector.h"
35 #include "fMatrix.h"
36 #include "lo-error.h"
37 #include "lo-lapack-proto.h"
38 #include "lo-qrupdate-proto.h"
39 #include "lu.h"
40 #include "oct-locbuf.h"
41 
42 namespace octave
43 {
44  namespace math
45  {
46  template <typename T>
47  lu<T>::lu (const T& l, const T& u,
48  const PermMatrix& p)
49  : a_fact (u), l_fact (l), ipvt (p.transpose ().col_perm_vec ())
50  {
51  if (l.columns () != u.rows ())
52  (*current_liboctave_error_handler) ("lu: dimension mismatch");
53  }
54 
55  template <typename T>
56  bool
57  lu<T>::packed (void) const
58  {
59  return l_fact.dims () == dim_vector ();
60  }
61 
62  template <typename T>
63  void
65  {
66  if (packed ())
67  {
68  l_fact = L ();
69  a_fact = U (); // FIXME: sub-optimal
70  ipvt = getp ();
71  }
72  }
73 
74  template <typename T>
75  T
76  lu<T>::L (void) const
77  {
78  if (packed ())
79  {
80  octave_idx_type a_nr = a_fact.rows ();
81  octave_idx_type a_nc = a_fact.cols ();
82  octave_idx_type mn = (a_nr < a_nc ? a_nr : a_nc);
83 
84  T l (a_nr, mn, ELT_T (0.0));
85 
86  for (octave_idx_type i = 0; i < a_nr; i++)
87  {
88  if (i < a_nc)
89  l.xelem (i, i) = 1.0;
90 
91  for (octave_idx_type j = 0; j < (i < a_nc ? i : a_nc); j++)
92  l.xelem (i, j) = a_fact.xelem (i, j);
93  }
94 
95  return l;
96  }
97  else
98  return l_fact;
99  }
100 
101  template <typename T>
102  T
103  lu<T>::U (void) const
104  {
105  if (packed ())
106  {
107  octave_idx_type a_nr = a_fact.rows ();
108  octave_idx_type a_nc = a_fact.cols ();
109  octave_idx_type mn = (a_nr < a_nc ? a_nr : a_nc);
110 
111  T u (mn, a_nc, ELT_T (0.0));
112 
113  for (octave_idx_type i = 0; i < mn; i++)
114  {
115  for (octave_idx_type j = i; j < a_nc; j++)
116  u.xelem (i, j) = a_fact.xelem (i, j);
117  }
118 
119  return u;
120  }
121  else
122  return a_fact;
123  }
124 
125  template <typename T>
126  T
127  lu<T>::Y (void) const
128  {
129  if (! packed ())
130  (*current_liboctave_error_handler)
131  ("lu: Y () not implemented for unpacked form");
132 
133  return a_fact;
134  }
135 
136  template <typename T>
138  lu<T>::getp (void) const
139  {
140  if (packed ())
141  {
142  octave_idx_type a_nr = a_fact.rows ();
143 
144  Array<octave_idx_type> pvt (dim_vector (a_nr, 1));
145 
146  for (octave_idx_type i = 0; i < a_nr; i++)
147  pvt.xelem (i) = i;
148 
149  for (octave_idx_type i = 0; i < ipvt.numel (); i++)
150  {
151  octave_idx_type k = ipvt.xelem (i);
152 
153  if (k != i)
154  {
155  octave_idx_type tmp = pvt.xelem (k);
156  pvt.xelem (k) = pvt.xelem (i);
157  pvt.xelem (i) = tmp;
158  }
159  }
160 
161  return pvt;
162  }
163  else
164  return ipvt;
165  }
166 
167  template <typename T>
168  PermMatrix
169  lu<T>::P (void) const
170  {
171  return PermMatrix (getp (), false);
172  }
173 
174  template <typename T>
176  lu<T>::P_vec (void) const
177  {
178  octave_idx_type a_nr = a_fact.rows ();
179 
180  ColumnVector p (a_nr);
181 
182  Array<octave_idx_type> pvt = getp ();
183 
184  for (octave_idx_type i = 0; i < a_nr; i++)
185  p.xelem (i) = static_cast<double> (pvt.xelem (i) + 1);
186 
187  return p;
188  }
189 
190  template <typename T>
191  bool
192  lu<T>::regular (void) const
193  {
194  bool retval = true;
195 
196  octave_idx_type k = std::min (a_fact.rows (), a_fact.columns ());
197 
198  for (octave_idx_type i = 0; i < k; i++)
199  {
200  if (a_fact(i, i) == ELT_T ())
201  {
202  retval = false;
203  break;
204  }
205  }
206 
207  return retval;
208  }
209 
210 #if ! defined (HAVE_QRUPDATE_LUU)
211 
212  template <typename T>
213  void
214  lu<T>::update (const VT&, const VT&)
215  {
216  (*current_liboctave_error_handler)
217  ("luupdate: support for qrupdate with LU updates "
218  "was unavailable or disabled when liboctave was built");
219  }
220 
221  template <typename T>
222  void
223  lu<T>::update (const T&, const T&)
224  {
225  (*current_liboctave_error_handler)
226  ("luupdate: support for qrupdate with LU updates "
227  "was unavailable or disabled when liboctave was built");
228  }
229 
230  template <typename T>
231  void
232  lu<T>::update_piv (const VT&, const VT&)
233  {
234  (*current_liboctave_error_handler)
235  ("luupdate: support for qrupdate with LU updates "
236  "was unavailable or disabled when liboctave was built");
237  }
238 
239  template <typename T>
240  void
241  lu<T>::update_piv (const T&, const T&)
242  {
243  (*current_liboctave_error_handler)
244  ("luupdate: support for qrupdate with LU updates "
245  "was unavailable or disabled when liboctave was built");
246  }
247 
248 #endif
249 
250  // Specializations.
251 
252  template <>
254  {
255  octave_idx_type a_nr = a.rows ();
256  octave_idx_type a_nc = a.cols ();
257  octave_idx_type mn = (a_nr < a_nc ? a_nr : a_nc);
258 
259  ipvt.resize (dim_vector (mn, 1));
260  octave_idx_type *pipvt = ipvt.fortran_vec ();
261 
262  a_fact = a;
263  double *tmp_data = a_fact.fortran_vec ();
264 
265  octave_idx_type info = 0;
266 
267  F77_XFCN (dgetrf, DGETRF, (a_nr, a_nc, tmp_data, a_nr, pipvt, info));
268 
269  for (octave_idx_type i = 0; i < mn; i++)
270  pipvt[i] -= 1;
271  }
272 
273 #if defined (HAVE_QRUPDATE_LUU)
274 
275  template <>
276  void
278  {
279  if (packed ())
280  unpack ();
281 
282  Matrix& l = l_fact;
283  Matrix& r = a_fact;
284 
285  octave_idx_type m = l.rows ();
286  octave_idx_type n = r.columns ();
287  octave_idx_type k = l.columns ();
288 
289  if (u.numel () != m || v.numel () != n)
290  (*current_liboctave_error_handler) ("luupdate: dimensions mismatch");
291 
292  ColumnVector utmp = u;
293  ColumnVector vtmp = v;
294  F77_XFCN (dlu1up, DLU1UP, (m, n, l.fortran_vec (), m, r.fortran_vec (), k,
295  utmp.fortran_vec (), vtmp.fortran_vec ()));
296  }
297 
298  template <>
299  void
300  lu<Matrix>::update (const Matrix& u, const Matrix& v)
301  {
302  if (packed ())
303  unpack ();
304 
305  Matrix& l = l_fact;
306  Matrix& r = a_fact;
307 
308  octave_idx_type m = l.rows ();
309  octave_idx_type n = r.columns ();
310  octave_idx_type k = l.columns ();
311 
312  if (u.rows () != m || v.rows () != n || u.cols () != v.cols ())
313  (*current_liboctave_error_handler) ("luupdate: dimensions mismatch");
314 
315  for (volatile octave_idx_type i = 0; i < u.cols (); i++)
316  {
317  ColumnVector utmp = u.column (i);
318  ColumnVector vtmp = v.column (i);
319  F77_XFCN (dlu1up, DLU1UP, (m, n, l.fortran_vec (),
320  m, r.fortran_vec (), k,
321  utmp.fortran_vec (), vtmp.fortran_vec ()));
322  }
323  }
324 
325  template <>
326  void
328  {
329  if (packed ())
330  unpack ();
331 
332  Matrix& l = l_fact;
333  Matrix& r = a_fact;
334 
335  octave_idx_type m = l.rows ();
336  octave_idx_type n = r.columns ();
337  octave_idx_type k = l.columns ();
338 
339  if (u.numel () != m || v.numel () != n)
340  (*current_liboctave_error_handler) ("luupdate: dimensions mismatch");
341 
342  ColumnVector utmp = u;
343  ColumnVector vtmp = v;
344  OCTAVE_LOCAL_BUFFER (double, w, m);
345  for (octave_idx_type i = 0; i < m; i++) ipvt(i) += 1; // increment
346  F77_XFCN (dlup1up, DLUP1UP, (m, n, l.fortran_vec (),
347  m, r.fortran_vec (), k,
348  ipvt.fortran_vec (),
349  utmp.data (), vtmp.data (), w));
350  for (octave_idx_type i = 0; i < m; i++) ipvt(i) -= 1; // decrement
351  }
352 
353  template <>
354  void
356  {
357  if (packed ())
358  unpack ();
359 
360  Matrix& l = l_fact;
361  Matrix& r = a_fact;
362 
363  octave_idx_type m = l.rows ();
364  octave_idx_type n = r.columns ();
365  octave_idx_type k = l.columns ();
366 
367  if (u.rows () != m || v.rows () != n || u.cols () != v.cols ())
368  (*current_liboctave_error_handler) ("luupdate: dimensions mismatch");
369 
370  OCTAVE_LOCAL_BUFFER (double, w, m);
371  for (octave_idx_type i = 0; i < m; i++) ipvt(i) += 1; // increment
372  for (volatile octave_idx_type i = 0; i < u.cols (); i++)
373  {
374  ColumnVector utmp = u.column (i);
375  ColumnVector vtmp = v.column (i);
376  F77_XFCN (dlup1up, DLUP1UP, (m, n, l.fortran_vec (),
377  m, r.fortran_vec (), k,
378  ipvt.fortran_vec (),
379  utmp.data (), vtmp.data (), w));
380  }
381  for (octave_idx_type i = 0; i < m; i++) ipvt(i) -= 1; // decrement
382  }
383 
384 #endif
385 
386  template <>
388  {
389  octave_idx_type a_nr = a.rows ();
390  octave_idx_type a_nc = a.cols ();
391  octave_idx_type mn = (a_nr < a_nc ? a_nr : a_nc);
392 
393  ipvt.resize (dim_vector (mn, 1));
394  octave_idx_type *pipvt = ipvt.fortran_vec ();
395 
396  a_fact = a;
397  float *tmp_data = a_fact.fortran_vec ();
398 
399  octave_idx_type info = 0;
400 
401  F77_XFCN (sgetrf, SGETRF, (a_nr, a_nc, tmp_data, a_nr, pipvt, info));
402 
403  for (octave_idx_type i = 0; i < mn; i++)
404  pipvt[i] -= 1;
405  }
406 
407 #if defined (HAVE_QRUPDATE_LUU)
408 
409  template <>
410  void
412  {
413  if (packed ())
414  unpack ();
415 
416  FloatMatrix& l = l_fact;
417  FloatMatrix& r = a_fact;
418 
419  octave_idx_type m = l.rows ();
420  octave_idx_type n = r.columns ();
421  octave_idx_type k = l.columns ();
422 
423  if (u.numel () != m || v.numel () != n)
424  (*current_liboctave_error_handler) ("luupdate: dimensions mismatch");
425 
426  FloatColumnVector utmp = u;
427  FloatColumnVector vtmp = v;
428  F77_XFCN (slu1up, SLU1UP, (m, n, l.fortran_vec (),
429  m, r.fortran_vec (), k,
430  utmp.fortran_vec (), vtmp.fortran_vec ()));
431  }
432 
433  template <>
434  void
436  {
437  if (packed ())
438  unpack ();
439 
440  FloatMatrix& l = l_fact;
441  FloatMatrix& r = a_fact;
442 
443  octave_idx_type m = l.rows ();
444  octave_idx_type n = r.columns ();
445  octave_idx_type k = l.columns ();
446 
447  if (u.rows () != m || v.rows () != n || u.cols () != v.cols ())
448  (*current_liboctave_error_handler) ("luupdate: dimensions mismatch");
449 
450  for (volatile octave_idx_type i = 0; i < u.cols (); i++)
451  {
452  FloatColumnVector utmp = u.column (i);
453  FloatColumnVector vtmp = v.column (i);
454  F77_XFCN (slu1up, SLU1UP, (m, n, l.fortran_vec (),
455  m, r.fortran_vec (), k,
456  utmp.fortran_vec (), vtmp.fortran_vec ()));
457  }
458  }
459 
460  template <>
461  void
463  const FloatColumnVector& v)
464  {
465  if (packed ())
466  unpack ();
467 
468  FloatMatrix& l = l_fact;
469  FloatMatrix& r = a_fact;
470 
471  octave_idx_type m = l.rows ();
472  octave_idx_type n = r.columns ();
473  octave_idx_type k = l.columns ();
474 
475  if (u.numel () != m || v.numel () != n)
476  (*current_liboctave_error_handler) ("luupdate: dimensions mismatch");
477 
478  FloatColumnVector utmp = u;
479  FloatColumnVector vtmp = v;
480  OCTAVE_LOCAL_BUFFER (float, w, m);
481  for (octave_idx_type i = 0; i < m; i++) ipvt(i) += 1; // increment
482  F77_XFCN (slup1up, SLUP1UP, (m, n, l.fortran_vec (),
483  m, r.fortran_vec (), k,
484  ipvt.fortran_vec (),
485  utmp.data (), vtmp.data (), w));
486  for (octave_idx_type i = 0; i < m; i++) ipvt(i) -= 1; // decrement
487  }
488 
489  template <>
490  void
492  {
493  if (packed ())
494  unpack ();
495 
496  FloatMatrix& l = l_fact;
497  FloatMatrix& r = a_fact;
498 
499  octave_idx_type m = l.rows ();
500  octave_idx_type n = r.columns ();
501  octave_idx_type k = l.columns ();
502 
503  if (u.rows () != m || v.rows () != n || u.cols () != v.cols ())
504  (*current_liboctave_error_handler) ("luupdate: dimensions mismatch");
505 
506  OCTAVE_LOCAL_BUFFER (float, w, m);
507  for (octave_idx_type i = 0; i < m; i++) ipvt(i) += 1; // increment
508  for (volatile octave_idx_type i = 0; i < u.cols (); i++)
509  {
510  FloatColumnVector utmp = u.column (i);
511  FloatColumnVector vtmp = v.column (i);
512  F77_XFCN (slup1up, SLUP1UP, (m, n, l.fortran_vec (),
513  m, r.fortran_vec (), k,
514  ipvt.fortran_vec (),
515  utmp.data (), vtmp.data (), w));
516  }
517  for (octave_idx_type i = 0; i < m; i++) ipvt(i) -= 1; // decrement
518  }
519 
520 #endif
521 
522  template <>
524  {
525  octave_idx_type a_nr = a.rows ();
526  octave_idx_type a_nc = a.cols ();
527  octave_idx_type mn = (a_nr < a_nc ? a_nr : a_nc);
528 
529  ipvt.resize (dim_vector (mn, 1));
530  octave_idx_type *pipvt = ipvt.fortran_vec ();
531 
532  a_fact = a;
533  Complex *tmp_data = a_fact.fortran_vec ();
534 
535  octave_idx_type info = 0;
536 
537  F77_XFCN (zgetrf, ZGETRF, (a_nr, a_nc, F77_DBLE_CMPLX_ARG (tmp_data), a_nr,
538  pipvt, info));
539 
540  for (octave_idx_type i = 0; i < mn; i++)
541  pipvt[i] -= 1;
542  }
543 
544 #if defined (HAVE_QRUPDATE_LUU)
545 
546  template <>
547  void
549  const ComplexColumnVector& v)
550  {
551  if (packed ())
552  unpack ();
553 
554  ComplexMatrix& l = l_fact;
555  ComplexMatrix& r = a_fact;
556 
557  octave_idx_type m = l.rows ();
558  octave_idx_type n = r.columns ();
559  octave_idx_type k = l.columns ();
560 
561  if (u.numel () != m || v.numel () != n)
562  (*current_liboctave_error_handler) ("luupdate: dimensions mismatch");
563 
564  ComplexColumnVector utmp = u;
565  ComplexColumnVector vtmp = v;
566  F77_XFCN (zlu1up, ZLU1UP, (m, n, F77_DBLE_CMPLX_ARG (l.fortran_vec ()), m,
569  F77_DBLE_CMPLX_ARG (vtmp.fortran_vec ())));
570  }
571 
572  template <>
573  void
575  {
576  if (packed ())
577  unpack ();
578 
579  ComplexMatrix& l = l_fact;
580  ComplexMatrix& r = a_fact;
581 
582  octave_idx_type m = l.rows ();
583  octave_idx_type n = r.columns ();
584  octave_idx_type k = l.columns ();
585 
586  if (u.rows () != m || v.rows () != n || u.cols () != v.cols ())
587  (*current_liboctave_error_handler) ("luupdate: dimensions mismatch");
588 
589  for (volatile octave_idx_type i = 0; i < u.cols (); i++)
590  {
591  ComplexColumnVector utmp = u.column (i);
592  ComplexColumnVector vtmp = v.column (i);
593  F77_XFCN (zlu1up, ZLU1UP, (m, n, F77_DBLE_CMPLX_ARG (l.fortran_vec ()),
594  m, F77_DBLE_CMPLX_ARG (r.fortran_vec ()), k,
596  F77_DBLE_CMPLX_ARG (vtmp.fortran_vec ())));
597  }
598  }
599 
600  template <>
601  void
603  const ComplexColumnVector& v)
604  {
605  if (packed ())
606  unpack ();
607 
608  ComplexMatrix& l = l_fact;
609  ComplexMatrix& r = a_fact;
610 
611  octave_idx_type m = l.rows ();
612  octave_idx_type n = r.columns ();
613  octave_idx_type k = l.columns ();
614 
615  if (u.numel () != m || v.numel () != n)
616  (*current_liboctave_error_handler) ("luupdate: dimensions mismatch");
617 
618  ComplexColumnVector utmp = u;
619  ComplexColumnVector vtmp = v;
621  for (octave_idx_type i = 0; i < m; i++) ipvt(i) += 1; // increment
622  F77_XFCN (zlup1up, ZLUP1UP, (m, n, F77_DBLE_CMPLX_ARG (l.fortran_vec ()),
623  m, F77_DBLE_CMPLX_ARG (r.fortran_vec ()), k,
624  ipvt.fortran_vec (),
625  F77_CONST_DBLE_CMPLX_ARG (utmp.data ()),
627  for (octave_idx_type i = 0; i < m; i++) ipvt(i) -= 1; // decrement
628  }
629 
630  template <>
631  void
633  {
634  if (packed ())
635  unpack ();
636 
637  ComplexMatrix& l = l_fact;
638  ComplexMatrix& r = a_fact;
639 
640  octave_idx_type m = l.rows ();
641  octave_idx_type n = r.columns ();
642  octave_idx_type k = l.columns ();
643 
644  if (u.rows () != m || v.rows () != n || u.cols () != v.cols ())
645  (*current_liboctave_error_handler) ("luupdate: dimensions mismatch");
646 
648  for (octave_idx_type i = 0; i < m; i++) ipvt(i) += 1; // increment
649  for (volatile octave_idx_type i = 0; i < u.cols (); i++)
650  {
651  ComplexColumnVector utmp = u.column (i);
652  ComplexColumnVector vtmp = v.column (i);
653  F77_XFCN (zlup1up, ZLUP1UP, (m, n, F77_DBLE_CMPLX_ARG (l.fortran_vec ()),
654  m, F77_DBLE_CMPLX_ARG (r.fortran_vec ()), k,
655  ipvt.fortran_vec (),
656  F77_CONST_DBLE_CMPLX_ARG (utmp.data ()),
658  }
659  for (octave_idx_type i = 0; i < m; i++) ipvt(i) -= 1; // decrement
660  }
661 
662 #endif
663 
664  template <>
666  {
667  octave_idx_type a_nr = a.rows ();
668  octave_idx_type a_nc = a.cols ();
669  octave_idx_type mn = (a_nr < a_nc ? a_nr : a_nc);
670 
671  ipvt.resize (dim_vector (mn, 1));
672  octave_idx_type *pipvt = ipvt.fortran_vec ();
673 
674  a_fact = a;
675  FloatComplex *tmp_data = a_fact.fortran_vec ();
676 
677  octave_idx_type info = 0;
678 
679  F77_XFCN (cgetrf, CGETRF, (a_nr, a_nc, F77_CMPLX_ARG (tmp_data), a_nr, pipvt,
680  info));
681 
682  for (octave_idx_type i = 0; i < mn; i++)
683  pipvt[i] -= 1;
684  }
685 
686 #if defined (HAVE_QRUPDATE_LUU)
687 
688  template <>
689  void
691  const FloatComplexColumnVector& v)
692  {
693  if (packed ())
694  unpack ();
695 
696  FloatComplexMatrix& l = l_fact;
697  FloatComplexMatrix& r = a_fact;
698 
699  octave_idx_type m = l.rows ();
700  octave_idx_type n = r.columns ();
701  octave_idx_type k = l.columns ();
702 
703  if (u.numel () == m && v.numel () == n)
704  {
706  FloatComplexColumnVector vtmp = v;
707  F77_XFCN (clu1up, CLU1UP, (m, n, F77_CMPLX_ARG (l.fortran_vec ()), m,
708  F77_CMPLX_ARG (r.fortran_vec ()), k,
709  F77_CMPLX_ARG (utmp.fortran_vec ()), F77_CMPLX_ARG (vtmp.fortran_vec ())));
710  }
711  else
712  (*current_liboctave_error_handler) ("luupdate: dimensions mismatch");
713  }
714 
715  template <>
716  void
718  const FloatComplexMatrix& v)
719  {
720  if (packed ())
721  unpack ();
722 
723  FloatComplexMatrix& l = l_fact;
724  FloatComplexMatrix& r = a_fact;
725 
726  octave_idx_type m = l.rows ();
727  octave_idx_type n = r.columns ();
728  octave_idx_type k = l.columns ();
729 
730  if (u.rows () != m || v.rows () != n || u.cols () != v.cols ())
731  (*current_liboctave_error_handler) ("luupdate: dimensions mismatch");
732 
733  for (volatile octave_idx_type i = 0; i < u.cols (); i++)
734  {
735  FloatComplexColumnVector utmp = u.column (i);
736  FloatComplexColumnVector vtmp = v.column (i);
737  F77_XFCN (clu1up, CLU1UP, (m, n, F77_CMPLX_ARG (l.fortran_vec ()),
738  m, F77_CMPLX_ARG (r.fortran_vec ()), k,
739  F77_CMPLX_ARG (utmp.fortran_vec ()), F77_CMPLX_ARG (vtmp.fortran_vec ())));
740  }
741  }
742 
743  template <>
744  void
746  const FloatComplexColumnVector& v)
747  {
748  if (packed ())
749  unpack ();
750 
751  FloatComplexMatrix& l = l_fact;
752  FloatComplexMatrix& r = a_fact;
753 
754  octave_idx_type m = l.rows ();
755  octave_idx_type n = r.columns ();
756  octave_idx_type k = l.columns ();
757 
758  if (u.numel () != m || v.numel () != n)
759  (*current_liboctave_error_handler) ("luupdate: dimensions mismatch");
760 
762  FloatComplexColumnVector vtmp = v;
764  for (octave_idx_type i = 0; i < m; i++) ipvt(i) += 1; // increment
765  F77_XFCN (clup1up, CLUP1UP, (m, n, F77_CMPLX_ARG (l.fortran_vec ()),
766  m, F77_CMPLX_ARG (r.fortran_vec ()), k,
767  ipvt.fortran_vec (),
768  F77_CONST_CMPLX_ARG (utmp.data ()), F77_CONST_CMPLX_ARG (vtmp.data ()),
769  F77_CMPLX_ARG (w)));
770  for (octave_idx_type i = 0; i < m; i++) ipvt(i) -= 1; // decrement
771  }
772 
773  template <>
774  void
776  const FloatComplexMatrix& v)
777  {
778  if (packed ())
779  unpack ();
780 
781  FloatComplexMatrix& l = l_fact;
782  FloatComplexMatrix& r = a_fact;
783 
784  octave_idx_type m = l.rows ();
785  octave_idx_type n = r.columns ();
786  octave_idx_type k = l.columns ();
787 
788  if (u.rows () != m || v.rows () != n || u.cols () != v.cols ())
789  (*current_liboctave_error_handler) ("luupdate: dimensions mismatch");
790 
792  for (octave_idx_type i = 0; i < m; i++) ipvt(i) += 1; // increment
793  for (volatile octave_idx_type i = 0; i < u.cols (); i++)
794  {
795  FloatComplexColumnVector utmp = u.column (i);
796  FloatComplexColumnVector vtmp = v.column (i);
797  F77_XFCN (clup1up, CLUP1UP, (m, n, F77_CMPLX_ARG (l.fortran_vec ()),
798  m, F77_CMPLX_ARG (r.fortran_vec ()), k,
799  ipvt.fortran_vec (),
800  F77_CONST_CMPLX_ARG (utmp.data ()), F77_CONST_CMPLX_ARG (vtmp.data ()),
801  F77_CMPLX_ARG (w)));
802  }
803  for (octave_idx_type i = 0; i < m; i++) ipvt(i) -= 1; // decrement
804  }
805 
806 #endif
807 
808  // Instantiations we need.
809 
810  template class lu<Matrix>;
811 
812  template class lu<FloatMatrix>;
813 
814  template class lu<ComplexMatrix>;
815 
816  template class lu<FloatComplexMatrix>;
817  }
818 }
bool packed(void) const
Definition: lu.cc:57
Octave interface to the compression and uncompression libraries.
Definition: aepbalance.cc:47
octave_idx_type numel(void) const
Number of elements in the array.
Definition: Array.h:363
OCTAVE_NORETURN liboctave_error_handler current_liboctave_error_handler
Definition: lo-error.c:38
#define F77_DBLE_CMPLX_ARG(x)
Definition: f77-fcn.h:345
for large enough k
Definition: lu.cc:606
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:382
u
Definition: lu.cc:138
#define F77_XFCN(f, F, args)
Definition: f77-fcn.h:52
octave_idx_type a_nc
Definition: sylvester.cc:72
Array< octave_idx_type > getp(void) const
Definition: lu.cc:138
octave_idx_type rows(void) const
Definition: Array.h:401
T L(void) const
Definition: lu.cc:76
calling an anonymous function involves an overhead quite comparable to the overhead of an m file function Passing a handle to a built in function is because the interpreter is not involved in the internal loop For a
Definition: cellfun.cc:398
bool regular(void) const
Definition: lu.cc:192
FloatComplexColumnVector column(octave_idx_type i) const
Definition: fCMatrix.cc:714
octave_idx_type a_nr
Definition: sylvester.cc:71
nd deftypefn *octave_map m
Definition: ov-struct.cc:2058
std::complex< double > w(std::complex< double > z, double relerr=0)
const T * data(void) const
Definition: Array.h:582
double tmp
Definition: data.cc:6300
PermMatrix P(void) const
Definition: lu.cc:169
octave_value retval
Definition: data.cc:6294
void update_piv(const VT &u, const VT &v)
Definition: dMatrix.h:37
FloatColumnVector column(octave_idx_type i) const
Definition: fMatrix.cc:429
T U(void) const
Definition: lu.cc:103
T & xelem(octave_idx_type n)
Definition: Array.h:455
lu(void)
Definition: lu.h:44
#define F77_CMPLX_ARG(x)
Definition: f77-fcn.h:339
#define F77_CONST_DBLE_CMPLX_ARG(x)
Definition: f77-fcn.h:348
T::element_type ELT_T
Definition: lu.h:42
=val(i)}if ode{val(i)}occurs in table i
Definition: lookup.cc:239
p
Definition: lu.cc:138
ColumnVector column(octave_idx_type i) const
Definition: dMatrix.cc:423
#define OCTAVE_LOCAL_BUFFER(T, buf, size)
Definition: oct-locbuf.h:200
void update(const VT &u, const VT &v)
#define F77_CONST_CMPLX_ARG(x)
Definition: f77-fcn.h:342
std::complex< float > FloatComplex
Definition: oct-cmplx.h:32
ColumnVector P_vec(void) const
Definition: lu.cc:176
void unpack(void)
Definition: lu.cc:64
std::complex< double > Complex
Definition: oct-cmplx.h:31
const T * fortran_vec(void) const
Definition: Array.h:584
octave_idx_type cols(void) const
Definition: Array.h:409
Vector representing the dimensions (size) of an Array.
Definition: dim-vector.h:87
T Y(void) const
Definition: lu.cc:127
Definition: mx-defs.h:75
ComplexColumnVector column(octave_idx_type i) const
Definition: CMatrix.cc:711
octave_idx_type columns(void) const
Definition: Array.h:410
charNDArray min(char d, const charNDArray &m)
Definition: chNDArray.cc:205