GNU Octave 7.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-2022 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
47namespace octave
48{
49 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 (void) const
69 {
70 return m_L.dims () == dim_vector ();
71 }
72
73 template <typename T>
74 void
75 lu<T>::unpack (void)
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 (void) const
97 {
98 if (packed ())
99 {
100 octave_idx_type a_nr = m_a_fact.rows ();
101 octave_idx_type a_nc = m_a_fact.columns ();
102 octave_idx_type mn = (a_nr < a_nc ? a_nr : a_nc);
103
104 T l (a_nr, mn, ELT_T (0.0));
105
106 for (octave_idx_type i = 0; i < a_nr; i++)
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 (void) const
124 {
125 if (packed ())
126 {
127 octave_idx_type a_nr = m_a_fact.rows ();
128 octave_idx_type a_nc = m_a_fact.columns ();
129 octave_idx_type mn = (a_nr < a_nc ? a_nr : a_nc);
130
131 T u (mn, a_nc, ELT_T (0.0));
132
133 for (octave_idx_type i = 0; i < mn; i++)
134 {
135 for (octave_idx_type j = i; j < a_nc; j++)
136 u.xelem (i, j) = m_a_fact.xelem (i, j);
137 }
138
139 return u;
140 }
141 else
142 return m_a_fact;
143 }
144
145 template <typename T>
146 T
147 lu<T>::Y (void) const
148 {
149 if (! packed ())
150 (*current_liboctave_error_handler)
151 ("lu: Y () not implemented for unpacked form");
152
153 return m_a_fact;
154 }
155
156 template <typename T>
157 Array<octave_idx_type>
158 lu<T>::getp (void) const
159 {
160 if (packed ())
161 {
162 octave_idx_type a_nr = m_a_fact.rows ();
163
164 Array<octave_idx_type> pvt (dim_vector (a_nr, 1));
165
166 for (octave_idx_type i = 0; i < a_nr; i++)
167 pvt.xelem (i) = i;
168
169 for (octave_idx_type i = 0; i < m_ipvt.numel (); i++)
170 {
171 octave_idx_type k = m_ipvt.xelem (i);
172
173 if (k != i)
174 {
175 octave_idx_type tmp = pvt.xelem (k);
176 pvt.xelem (k) = pvt.xelem (i);
177 pvt.xelem (i) = tmp;
178 }
179 }
180
181 return pvt;
182 }
183 else
184 return m_ipvt;
185 }
186
187 template <typename T>
188 PermMatrix
189 lu<T>::P (void) const
190 {
191 return PermMatrix (getp (), false);
192 }
193
194 template <typename T>
195 ColumnVector
196 lu<T>::P_vec (void) const
197 {
198 octave_idx_type a_nr = m_a_fact.rows ();
199
200 ColumnVector p (a_nr);
201
202 Array<octave_idx_type> pvt = getp ();
203
204 for (octave_idx_type i = 0; i < a_nr; i++)
205 p.xelem (i) = static_cast<double> (pvt.xelem (i) + 1);
206
207 return p;
208 }
209
210 template <typename T>
211 bool
212 lu<T>::regular (void) const
213 {
214 bool retval = true;
215
216 octave_idx_type k = std::min (m_a_fact.rows (), m_a_fact.columns ());
217
218 for (octave_idx_type i = 0; i < k; i++)
219 {
220 if (m_a_fact(i, i) == ELT_T ())
221 {
222 retval = false;
223 break;
224 }
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}
min
charNDArray min(char d, const charNDArray &m)
Definition: chNDArray.cc:207
Array::xelem
T & xelem(octave_idx_type n)
Size of the specified dimension.
Definition: Array.h:504
Array::numel
octave_idx_type numel(void) const
Number of elements in the array.
Definition: Array.h:411
Array::rows
octave_idx_type rows(void) const
Definition: Array.h:449
Array::columns
octave_idx_type columns(void) const
Definition: Array.h:458
Array::data
const T * data(void) const
Size of the specified dimension.
Definition: Array.h:616
Array::fortran_vec
OCTARRAY_API T * fortran_vec(void)
Size of the specified dimension.
Definition: Array.cc:1744
ComplexMatrix::column
OCTAVE_API ComplexColumnVector column(octave_idx_type i) const
Definition: CMatrix.cc:708
FloatComplexMatrix::column
OCTAVE_API FloatComplexColumnVector column(octave_idx_type i) const
Definition: fCMatrix.cc:711
FloatMatrix::column
OCTAVE_API FloatColumnVector column(octave_idx_type i) const
Definition: fMatrix.cc:428
Matrix
Definition: dMatrix.h:42
Matrix::column
OCTAVE_API 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: mx-defs.h:51
octave::math::lu::update
OCTAVE_API void update(const VT &u, const VT &v)
octave::math::lu::lu
lu(void)
Definition: lu.h:48
octave::math::lu::U
OCTAVE_API T U(void) const
Definition: lu.cc:123
octave::math::lu::regular
OCTAVE_API bool regular(void) const
Definition: lu.cc:212
octave::math::lu::P
OCTAVE_API PermMatrix P(void) const
Definition: lu.cc:189
octave::math::lu::Y
OCTAVE_API T Y(void) const
Definition: lu.cc:147
octave::math::lu::update_piv
OCTAVE_API void update_piv(const VT &u, const VT &v)
octave::math::lu::unpack
OCTAVE_API void unpack(void)
Definition: lu.cc:75
octave::math::lu::packed
OCTAVE_API bool packed(void) const
Definition: lu.cc:68
octave::math::lu::ELT_T
T::element_type ELT_T
Definition: lu.h:46
octave::math::lu::P_vec
OCTAVE_API ColumnVector P_vec(void) const
Definition: lu.cc:196
octave::math::lu::getp
OCTAVE_API Array< octave_idx_type > getp(void) const
Definition: lu.cc:158
octave::math::lu::L
OCTAVE_API 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: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.in.cc:55
PermMatrix
class OCTAVE_API PermMatrix
Definition: mx-fwd.h:64
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
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:391