GNU Octave 11.1.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 
Loading...
Searching...
No Matches
lu.cc
Go to the documentation of this file.
1////////////////////////////////////////////////////////////////////////
2//
3// Copyright (C) 1996-2026 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 "lapack-proto.h"
42#include "lu.h"
43#include "oct-error.h"
44#include "oct-locbuf.h"
45#include "qrupdate-proto.h"
46
47OCTAVE_BEGIN_NAMESPACE(octave)
48OCTAVE_BEGIN_NAMESPACE(math)
49
50// FIXME: PermMatrix::col_perm_vec returns Array<octave_idx_type>
51// but m_ipvt is an Array<octave_f77_int_type>. This could cause
52// trouble for large arrays if octave_f77_int_type is 32-bits but
53// octave_idx_type is 64. Since this constructor is called from
54// Fluupdate, it could be given values that are out of range. We
55// should ensure that the values are within range here.
56
57template <typename T>
58lu<T>::lu (const T& l, const T& u, const PermMatrix& p)
59 : m_a_fact (u), m_L (l), m_ipvt (p.transpose ().col_perm_vec ())
60{
61 if (l.columns () != u.rows ())
62 (*current_liboctave_error_handler) ("lu: dimension mismatch");
63}
64
65template <typename T>
66bool
67lu<T>::packed () const
68{
69 return m_L.dims () == dim_vector ();
70}
71
72template <typename T>
73void
74lu<T>::unpack ()
75{
76 if (packed ())
77 {
78 m_L = L ();
79 m_a_fact = U (); // FIXME: sub-optimal
80
81 // FIXME: getp returns Array<octave_idx_type> but m_ipvt is
82 // Array<octave_f77_int_type>. However, getp produces its
83 // result from a valid m_ipvt array so validation should not be
84 // necessary. OTOH, it might be better to have a version of
85 // getp that doesn't cause us to convert from
86 // Array<octave_f77_int_type> to Array<octave_idx_type> and
87 // back again.
88
89 m_ipvt = getp ();
90 }
91}
92
93template <typename T>
94T
95lu<T>::L () const
96{
97 if (packed ())
98 {
99 octave_idx_type a_nr = m_a_fact.rows ();
100 octave_idx_type a_nc = m_a_fact.columns ();
101 octave_idx_type mn = (a_nr < a_nc ? a_nr : a_nc);
102
103 T l (a_nr, mn, ELT_T (0.0));
104
105 for (octave_idx_type i = 0; i < a_nr; i++)
106 {
107 if (i < a_nc)
108 l.xelem (i, i) = 1.0;
109
110 for (octave_idx_type j = 0; j < (i < a_nc ? i : a_nc); j++)
111 l.xelem (i, j) = m_a_fact.xelem (i, j);
112 }
113
114 return l;
115 }
116 else
117 return m_L;
118}
119
120template <typename T>
121T
122lu<T>::U () const
123{
124 if (packed ())
125 {
126 octave_idx_type a_nr = m_a_fact.rows ();
127 octave_idx_type a_nc = m_a_fact.columns ();
128 octave_idx_type mn = (a_nr < a_nc ? a_nr : a_nc);
129
130 T u (mn, a_nc, ELT_T (0.0));
131
132 for (octave_idx_type i = 0; i < mn; i++)
133 {
134 for (octave_idx_type j = i; j < a_nc; j++)
135 u.xelem (i, j) = m_a_fact.xelem (i, j);
136 }
137
138 return u;
139 }
140 else
141 return m_a_fact;
142}
143
144template <typename T>
145T
146lu<T>::Y () const
147{
148 if (! packed ())
149 (*current_liboctave_error_handler)
150 ("lu: Y () not implemented for unpacked form");
151
152 return m_a_fact;
153}
154
155template <typename T>
156Array<octave_idx_type>
157lu<T>::getp () const
158{
159 if (packed ())
160 {
161 octave_idx_type a_nr = m_a_fact.rows ();
162
163 Array<octave_idx_type> pvt (dim_vector (a_nr, 1));
164
165 for (octave_idx_type i = 0; i < a_nr; i++)
166 pvt.xelem (i) = i;
167
168 for (octave_idx_type i = 0; i < m_ipvt.numel (); i++)
169 {
170 octave_idx_type k = m_ipvt.xelem (i);
171
172 if (k != i)
173 {
174 octave_idx_type tmp = pvt.xelem (k);
175 pvt.xelem (k) = pvt.xelem (i);
176 pvt.xelem (i) = tmp;
177 }
178 }
179
180 return pvt;
181 }
182 else
183 return m_ipvt;
184}
185
186template <typename T>
187PermMatrix
188lu<T>::P () const
189{
190 return PermMatrix (getp (), false);
191}
192
193template <typename T>
194ColumnVector
195lu<T>::P_vec () const
196{
197 octave_idx_type a_nr = m_a_fact.rows ();
198
199 ColumnVector p (a_nr);
200
201 Array<octave_idx_type> pvt = getp ();
202
203 for (octave_idx_type i = 0; i < a_nr; i++)
204 p.xelem (i) = static_cast<double> (pvt.xelem (i) + 1);
205
206 return p;
207}
208
209template <typename T>
210bool
211lu<T>::regular () const
212{
213 bool retval = true;
214
215 octave_idx_type k = std::min (m_a_fact.rows (), m_a_fact.columns ());
216
217 for (octave_idx_type i = 0; i < k; i++)
218 {
219 if (m_a_fact(i, i) == ELT_T ())
220 {
221 retval = false;
222 break;
223 }
224 }
225
226 return retval;
227}
228
229#if ! defined (HAVE_QRUPDATE_LUU)
230
231template <typename T>
232void
233lu<T>::update (const VT&, const VT&)
234{
235 (*current_liboctave_error_handler)
236 ("luupdate: support for qrupdate with LU updates "
237 "was unavailable or disabled when liboctave was built");
238}
239
240template <typename T>
241void
242lu<T>::update (const T&, const T&)
243{
244 (*current_liboctave_error_handler)
245 ("luupdate: support for qrupdate with LU updates "
246 "was unavailable or disabled when liboctave was built");
247}
248
249template <typename T>
250void
251lu<T>::update_piv (const VT&, const VT&)
252{
253 (*current_liboctave_error_handler)
254 ("luupdate: support for qrupdate with LU updates "
255 "was unavailable or disabled when liboctave was built");
256}
257
258template <typename T>
259void
260lu<T>::update_piv (const T&, const T&)
261{
262 (*current_liboctave_error_handler)
263 ("luupdate: support for qrupdate with LU updates "
264 "was unavailable or disabled when liboctave was built");
265}
266
267#endif
268
269// Specializations.
270
271template <>
272OCTAVE_API
273lu<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 m_ipvt.resize (dim_vector (mn, 1));
280 F77_INT *pipvt = m_ipvt.rwdata ();
281
282 m_a_fact = a;
283 double *tmp_data = m_a_fact.rwdata ();
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
295template <>
296OCTAVE_API void
297lu<Matrix>::update (const ColumnVector& u, const ColumnVector& v)
298{
299 if (packed ())
300 unpack ();
301
302 Matrix& l = m_L;
303 Matrix& r = m_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.rwdata (), m, r.rwdata (),
318 k, utmp.rwdata (), vtmp.rwdata ()));
319}
320
321template <>
322OCTAVE_API void
323lu<Matrix>::update (const Matrix& u, const Matrix& v)
324{
325 if (packed ())
326 unpack ();
327
328 Matrix& l = m_L;
329 Matrix& r = m_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 (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.rwdata (),
349 m, r.rwdata (), k,
350 utmp.rwdata (), vtmp.rwdata ()));
351 }
352}
353
354template <>
355OCTAVE_API void
356lu<Matrix>::update_piv (const ColumnVector& u, const ColumnVector& v)
357{
358 if (packed ())
359 unpack ();
360
361 Matrix& l = m_L;
362 Matrix& r = m_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++) m_ipvt(i) += 1; // increment
378 F77_XFCN (dlup1up, DLUP1UP, (m, n, l.rwdata (),
379 m, r.rwdata (), k,
380 m_ipvt.rwdata (),
381 utmp.data (), vtmp.data (), w));
382 for (F77_INT i = 0; i < m; i++) m_ipvt(i) -= 1; // decrement
383}
384
385template <>
386OCTAVE_API void
387lu<Matrix>::update_piv (const Matrix& u, const Matrix& v)
388{
389 if (packed ())
390 unpack ();
391
392 Matrix& l = m_L;
393 Matrix& r = m_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++) m_ipvt(i) += 1; // increment
410 for (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.rwdata (),
415 m, r.rwdata (), k,
416 m_ipvt.rwdata (),
417 utmp.data (), vtmp.data (), w));
418 }
419 for (F77_INT i = 0; i < m; i++) m_ipvt(i) -= 1; // decrement
420}
421
422#endif
423
424template <>
425OCTAVE_API
426lu<FloatMatrix>::lu (const FloatMatrix& a)
427{
428 F77_INT a_nr = to_f77_int (a.rows ());
429 F77_INT a_nc = to_f77_int (a.columns ());
430 F77_INT mn = (a_nr < a_nc ? a_nr : a_nc);
431
432 m_ipvt.resize (dim_vector (mn, 1));
433 F77_INT *pipvt = m_ipvt.rwdata ();
434
435 m_a_fact = a;
436 float *tmp_data = m_a_fact.rwdata ();
437
438 F77_INT info = 0;
439
440 F77_XFCN (sgetrf, SGETRF, (a_nr, a_nc, tmp_data, a_nr, pipvt, info));
441
442 for (F77_INT i = 0; i < mn; i++)
443 pipvt[i] -= 1;
444}
445
446#if defined (HAVE_QRUPDATE_LUU)
447
448template <>
449OCTAVE_API void
450lu<FloatMatrix>::update (const FloatColumnVector& u,
451 const FloatColumnVector& v)
452{
453 if (packed ())
454 unpack ();
455
456 FloatMatrix& l = m_L;
457 FloatMatrix& r = m_a_fact;
458
459 F77_INT m = to_f77_int (l.rows ());
460 F77_INT n = to_f77_int (r.columns ());
461 F77_INT k = to_f77_int (l.columns ());
462
463 F77_INT u_nel = to_f77_int (u.numel ());
464 F77_INT v_nel = to_f77_int (v.numel ());
465
466 if (u_nel != m || v_nel != n)
467 (*current_liboctave_error_handler) ("luupdate: dimensions mismatch");
468
469 FloatColumnVector utmp = u;
470 FloatColumnVector vtmp = v;
471 F77_XFCN (slu1up, SLU1UP, (m, n, l.rwdata (),
472 m, r.rwdata (), k,
473 utmp.rwdata (), vtmp.rwdata ()));
474}
475
476template <>
477OCTAVE_API void
478lu<FloatMatrix>::update (const FloatMatrix& u, const FloatMatrix& v)
479{
480 if (packed ())
481 unpack ();
482
483 FloatMatrix& l = m_L;
484 FloatMatrix& r = m_a_fact;
485
486 F77_INT m = to_f77_int (l.rows ());
487 F77_INT n = to_f77_int (r.columns ());
488 F77_INT k = to_f77_int (l.columns ());
489
490 F77_INT u_nr = to_f77_int (u.rows ());
491 F77_INT u_nc = to_f77_int (u.columns ());
492
493 F77_INT v_nr = to_f77_int (v.rows ());
494 F77_INT v_nc = to_f77_int (v.columns ());
495
496 if (u_nr != m || v_nr != n || u_nc != v_nc)
497 (*current_liboctave_error_handler) ("luupdate: dimensions mismatch");
498
499 for (F77_INT i = 0; i < u_nc; i++)
500 {
501 FloatColumnVector utmp = u.column (i);
502 FloatColumnVector vtmp = v.column (i);
503 F77_XFCN (slu1up, SLU1UP, (m, n, l.rwdata (),
504 m, r.rwdata (), k,
505 utmp.rwdata (), vtmp.rwdata ()));
506 }
507}
508
509template <>
510OCTAVE_API void
511lu<FloatMatrix>::update_piv (const FloatColumnVector& u,
512 const FloatColumnVector& v)
513{
514 if (packed ())
515 unpack ();
516
517 FloatMatrix& l = m_L;
518 FloatMatrix& r = m_a_fact;
519
520 F77_INT m = to_f77_int (l.rows ());
521 F77_INT n = to_f77_int (r.columns ());
522 F77_INT k = to_f77_int (l.columns ());
523
524 F77_INT u_nel = to_f77_int (u.numel ());
525 F77_INT v_nel = to_f77_int (v.numel ());
526
527 if (u_nel != m || v_nel != n)
528 (*current_liboctave_error_handler) ("luupdate: dimensions mismatch");
529
530 FloatColumnVector utmp = u;
531 FloatColumnVector vtmp = v;
532 OCTAVE_LOCAL_BUFFER (float, w, m);
533 for (F77_INT i = 0; i < m; i++) m_ipvt(i) += 1; // increment
534 F77_XFCN (slup1up, SLUP1UP, (m, n, l.rwdata (),
535 m, r.rwdata (), k,
536 m_ipvt.rwdata (),
537 utmp.data (), vtmp.data (), w));
538 for (F77_INT i = 0; i < m; i++) m_ipvt(i) -= 1; // decrement
539}
540
541template <>
542OCTAVE_API void
543lu<FloatMatrix>::update_piv (const FloatMatrix& u, const FloatMatrix& v)
544{
545 if (packed ())
546 unpack ();
547
548 FloatMatrix& l = m_L;
549 FloatMatrix& r = m_a_fact;
550
551 F77_INT m = to_f77_int (l.rows ());
552 F77_INT n = to_f77_int (r.columns ());
553 F77_INT k = to_f77_int (l.columns ());
554
555 F77_INT u_nr = to_f77_int (u.rows ());
556 F77_INT u_nc = to_f77_int (u.columns ());
557
558 F77_INT v_nr = to_f77_int (v.rows ());
559 F77_INT v_nc = to_f77_int (v.columns ());
560
561 if (u_nr != m || v_nr != n || u_nc != v_nc)
562 (*current_liboctave_error_handler) ("luupdate: dimensions mismatch");
563
564 OCTAVE_LOCAL_BUFFER (float, w, m);
565 for (F77_INT i = 0; i < m; i++) m_ipvt(i) += 1; // increment
566 for (F77_INT i = 0; i < u_nc; i++)
567 {
568 FloatColumnVector utmp = u.column (i);
569 FloatColumnVector vtmp = v.column (i);
570 F77_XFCN (slup1up, SLUP1UP, (m, n, l.rwdata (),
571 m, r.rwdata (), k,
572 m_ipvt.rwdata (),
573 utmp.data (), vtmp.data (), w));
574 }
575 for (F77_INT i = 0; i < m; i++) m_ipvt(i) -= 1; // decrement
576}
577
578#endif
579
580template <>
581OCTAVE_API
582lu<ComplexMatrix>::lu (const ComplexMatrix& a)
583{
584 F77_INT a_nr = to_f77_int (a.rows ());
585 F77_INT a_nc = to_f77_int (a.columns ());
586 F77_INT mn = (a_nr < a_nc ? a_nr : a_nc);
587
588 m_ipvt.resize (dim_vector (mn, 1));
589 F77_INT *pipvt = m_ipvt.rwdata ();
590
591 m_a_fact = a;
592 Complex *tmp_data = m_a_fact.rwdata ();
593
594 F77_INT info = 0;
595
596 F77_XFCN (zgetrf, ZGETRF, (a_nr, a_nc, F77_DBLE_CMPLX_ARG (tmp_data),
597 a_nr, pipvt, info));
598
599 for (F77_INT i = 0; i < mn; i++)
600 pipvt[i] -= 1;
601}
602
603#if defined (HAVE_QRUPDATE_LUU)
604
605template <>
606OCTAVE_API void
607lu<ComplexMatrix>::update (const ComplexColumnVector& u,
608 const ComplexColumnVector& v)
609{
610 if (packed ())
611 unpack ();
612
613 ComplexMatrix& l = m_L;
614 ComplexMatrix& r = m_a_fact;
615
616 F77_INT m = to_f77_int (l.rows ());
617 F77_INT n = to_f77_int (r.columns ());
618 F77_INT k = to_f77_int (l.columns ());
619
620 F77_INT u_nel = to_f77_int (u.numel ());
621 F77_INT v_nel = to_f77_int (v.numel ());
622
623 if (u_nel != m || v_nel != n)
624 (*current_liboctave_error_handler) ("luupdate: dimensions mismatch");
625
626 ComplexColumnVector utmp = u;
627 ComplexColumnVector vtmp = v;
628 F77_XFCN (zlu1up, ZLU1UP, (m, n, F77_DBLE_CMPLX_ARG (l.rwdata ()), m,
629 F77_DBLE_CMPLX_ARG (r.rwdata ()), k,
630 F77_DBLE_CMPLX_ARG (utmp.rwdata ()),
631 F77_DBLE_CMPLX_ARG (vtmp.rwdata ())));
632}
633
634template <>
635OCTAVE_API void
636lu<ComplexMatrix>::update (const ComplexMatrix& u, const ComplexMatrix& v)
637{
638 if (packed ())
639 unpack ();
640
641 ComplexMatrix& l = m_L;
642 ComplexMatrix& r = m_a_fact;
643
644 F77_INT m = to_f77_int (l.rows ());
645 F77_INT n = to_f77_int (r.columns ());
646 F77_INT k = to_f77_int (l.columns ());
647
648 F77_INT u_nr = to_f77_int (u.rows ());
649 F77_INT u_nc = to_f77_int (u.columns ());
650
651 F77_INT v_nr = to_f77_int (v.rows ());
652 F77_INT v_nc = to_f77_int (v.columns ());
653
654 if (u_nr != m || v_nr != n || u_nc != v_nc)
655 (*current_liboctave_error_handler) ("luupdate: dimensions mismatch");
656
657 for (F77_INT i = 0; i < u_nc; i++)
658 {
659 ComplexColumnVector utmp = u.column (i);
660 ComplexColumnVector vtmp = v.column (i);
661 F77_XFCN (zlu1up, ZLU1UP, (m, n,
662 F77_DBLE_CMPLX_ARG (l.rwdata ()),
663 m, F77_DBLE_CMPLX_ARG (r.rwdata ()),
664 k,
665 F77_DBLE_CMPLX_ARG (utmp.rwdata ()),
666 F77_DBLE_CMPLX_ARG (vtmp.rwdata ())));
667 }
668}
669
670template <>
671OCTAVE_API void
672lu<ComplexMatrix>::update_piv (const ComplexColumnVector& u,
673 const ComplexColumnVector& v)
674{
675 if (packed ())
676 unpack ();
677
678 ComplexMatrix& l = m_L;
679 ComplexMatrix& r = m_a_fact;
680
681 F77_INT m = to_f77_int (l.rows ());
682 F77_INT n = to_f77_int (r.columns ());
683 F77_INT k = to_f77_int (l.columns ());
684
685 F77_INT u_nel = to_f77_int (u.numel ());
686 F77_INT v_nel = to_f77_int (v.numel ());
687
688 if (u_nel != m || v_nel != n)
689 (*current_liboctave_error_handler) ("luupdate: dimensions mismatch");
690
691 ComplexColumnVector utmp = u;
692 ComplexColumnVector vtmp = v;
693 OCTAVE_LOCAL_BUFFER (Complex, w, m);
694 for (F77_INT i = 0; i < m; i++) m_ipvt(i) += 1; // increment
695 F77_XFCN (zlup1up, ZLUP1UP, (m, n, F77_DBLE_CMPLX_ARG (l.rwdata ()),
696 m, F77_DBLE_CMPLX_ARG (r.rwdata ()), k,
697 m_ipvt.rwdata (),
698 F77_CONST_DBLE_CMPLX_ARG (utmp.data ()),
699 F77_CONST_DBLE_CMPLX_ARG (vtmp.data ()),
700 F77_DBLE_CMPLX_ARG (w)));
701 for (F77_INT i = 0; i < m; i++) m_ipvt(i) -= 1; // decrement
702}
703
704template <>
705OCTAVE_API void
706lu<ComplexMatrix>::update_piv (const ComplexMatrix& u,
707 const ComplexMatrix& v)
708{
709 if (packed ())
710 unpack ();
711
712 ComplexMatrix& l = m_L;
713 ComplexMatrix& r = m_a_fact;
714
715 F77_INT m = to_f77_int (l.rows ());
716 F77_INT n = to_f77_int (r.columns ());
717 F77_INT k = to_f77_int (l.columns ());
718
719 F77_INT u_nr = to_f77_int (u.rows ());
720 F77_INT u_nc = to_f77_int (u.columns ());
721
722 F77_INT v_nr = to_f77_int (v.rows ());
723 F77_INT v_nc = to_f77_int (v.columns ());
724
725 if (u_nr != m || v_nr != n || u_nc != v_nc)
726 (*current_liboctave_error_handler) ("luupdate: dimensions mismatch");
727
728 OCTAVE_LOCAL_BUFFER (Complex, w, m);
729 for (F77_INT i = 0; i < m; i++) m_ipvt(i) += 1; // increment
730 for (F77_INT i = 0; i < u_nc; i++)
731 {
732 ComplexColumnVector utmp = u.column (i);
733 ComplexColumnVector vtmp = v.column (i);
734 F77_XFCN (zlup1up, ZLUP1UP, (m, n,
735 F77_DBLE_CMPLX_ARG (l.rwdata ()),
736 m,
737 F77_DBLE_CMPLX_ARG (r.rwdata ()),
738 k, m_ipvt.rwdata (),
739 F77_CONST_DBLE_CMPLX_ARG (utmp.data ()),
740 F77_CONST_DBLE_CMPLX_ARG (vtmp.data ()),
741 F77_DBLE_CMPLX_ARG (w)));
742 }
743 for (F77_INT i = 0; i < m; i++) m_ipvt(i) -= 1; // decrement
744}
745
746#endif
747
748template <>
749OCTAVE_API
750lu<FloatComplexMatrix>::lu (const FloatComplexMatrix& a)
751{
752 F77_INT a_nr = to_f77_int (a.rows ());
753 F77_INT a_nc = to_f77_int (a.columns ());
754 F77_INT mn = (a_nr < a_nc ? a_nr : a_nc);
755
756 m_ipvt.resize (dim_vector (mn, 1));
757 F77_INT *pipvt = m_ipvt.rwdata ();
758
759 m_a_fact = a;
760 FloatComplex *tmp_data = m_a_fact.rwdata ();
761
762 F77_INT info = 0;
763
764 F77_XFCN (cgetrf, CGETRF, (a_nr, a_nc, F77_CMPLX_ARG (tmp_data), a_nr,
765 pipvt, info));
766
767 for (F77_INT i = 0; i < mn; i++)
768 pipvt[i] -= 1;
769}
770
771#if defined (HAVE_QRUPDATE_LUU)
772
773template <>
774OCTAVE_API void
775lu<FloatComplexMatrix>::update (const FloatComplexColumnVector& u,
776 const FloatComplexColumnVector& v)
777{
778 if (packed ())
779 unpack ();
780
781 FloatComplexMatrix& l = m_L;
782 FloatComplexMatrix& r = m_a_fact;
783
784 F77_INT m = to_f77_int (l.rows ());
785 F77_INT n = to_f77_int (r.columns ());
786 F77_INT k = to_f77_int (l.columns ());
787
788 F77_INT u_nel = to_f77_int (u.numel ());
789 F77_INT v_nel = to_f77_int (v.numel ());
790
791 if (u_nel != m || v_nel != n)
792 (*current_liboctave_error_handler) ("luupdate: dimensions mismatch");
793
794 FloatComplexColumnVector utmp = u;
795 FloatComplexColumnVector vtmp = v;
796 F77_XFCN (clu1up, CLU1UP, (m, n, F77_CMPLX_ARG (l.rwdata ()), m,
797 F77_CMPLX_ARG (r.rwdata ()), k,
798 F77_CMPLX_ARG (utmp.rwdata ()),
799 F77_CMPLX_ARG (vtmp.rwdata ())));
800}
801
802template <>
803OCTAVE_API void
804lu<FloatComplexMatrix>::update (const FloatComplexMatrix& u,
805 const FloatComplexMatrix& v)
806{
807 if (packed ())
808 unpack ();
809
810 FloatComplexMatrix& l = m_L;
811 FloatComplexMatrix& r = m_a_fact;
812
813 F77_INT m = to_f77_int (l.rows ());
814 F77_INT n = to_f77_int (r.columns ());
815 F77_INT k = to_f77_int (l.columns ());
816
817 F77_INT u_nr = to_f77_int (u.rows ());
818 F77_INT u_nc = to_f77_int (u.columns ());
819
820 F77_INT v_nr = to_f77_int (v.rows ());
821 F77_INT v_nc = to_f77_int (v.columns ());
822
823 if (u_nr != m || v_nr != n || u_nc != v_nc)
824 (*current_liboctave_error_handler) ("luupdate: dimensions mismatch");
825
826 for (F77_INT i = 0; i < u_nc; i++)
827 {
828 FloatComplexColumnVector utmp = u.column (i);
829 FloatComplexColumnVector vtmp = v.column (i);
830 F77_XFCN (clu1up, CLU1UP, (m, n, F77_CMPLX_ARG (l.rwdata ()),
831 m, F77_CMPLX_ARG (r.rwdata ()), k,
832 F77_CMPLX_ARG (utmp.rwdata ()),
833 F77_CMPLX_ARG (vtmp.rwdata ())));
834 }
835}
836
837template <>
838OCTAVE_API void
839lu<FloatComplexMatrix>::update_piv (const FloatComplexColumnVector& u,
840 const FloatComplexColumnVector& v)
841{
842 if (packed ())
843 unpack ();
844
845 FloatComplexMatrix& l = m_L;
846 FloatComplexMatrix& r = m_a_fact;
847
848 F77_INT m = to_f77_int (l.rows ());
849 F77_INT n = to_f77_int (r.columns ());
850 F77_INT k = to_f77_int (l.columns ());
851
852 F77_INT u_nel = to_f77_int (u.numel ());
853 F77_INT v_nel = to_f77_int (v.numel ());
854
855 if (u_nel != m || v_nel != n)
856 (*current_liboctave_error_handler) ("luupdate: dimensions mismatch");
857
858 FloatComplexColumnVector utmp = u;
859 FloatComplexColumnVector vtmp = v;
860 OCTAVE_LOCAL_BUFFER (FloatComplex, w, m);
861 for (F77_INT i = 0; i < m; i++) m_ipvt(i) += 1; // increment
862 F77_XFCN (clup1up, CLUP1UP, (m, n, F77_CMPLX_ARG (l.rwdata ()),
863 m, F77_CMPLX_ARG (r.rwdata ()), k,
864 m_ipvt.rwdata (),
865 F77_CONST_CMPLX_ARG (utmp.data ()),
866 F77_CONST_CMPLX_ARG (vtmp.data ()),
867 F77_CMPLX_ARG (w)));
868 for (F77_INT i = 0; i < m; i++) m_ipvt(i) -= 1; // decrement
869}
870
871template <>
872OCTAVE_API void
873lu<FloatComplexMatrix>::update_piv (const FloatComplexMatrix& u,
874 const FloatComplexMatrix& v)
875{
876 if (packed ())
877 unpack ();
878
879 FloatComplexMatrix& l = m_L;
880 FloatComplexMatrix& r = m_a_fact;
881
882 F77_INT m = to_f77_int (l.rows ());
883 F77_INT n = to_f77_int (r.columns ());
884 F77_INT k = to_f77_int (l.columns ());
885
886 F77_INT u_nr = to_f77_int (u.rows ());
887 F77_INT u_nc = to_f77_int (u.columns ());
888
889 F77_INT v_nr = to_f77_int (v.rows ());
890 F77_INT v_nc = to_f77_int (v.columns ());
891
892 if (u_nr != m || v_nr != n || u_nc != v_nc)
893 (*current_liboctave_error_handler) ("luupdate: dimensions mismatch");
894
895 OCTAVE_LOCAL_BUFFER (FloatComplex, w, m);
896 for (F77_INT i = 0; i < m; i++) m_ipvt(i) += 1; // increment
897 for (F77_INT i = 0; i < u_nc; i++)
898 {
899 FloatComplexColumnVector utmp = u.column (i);
900 FloatComplexColumnVector vtmp = v.column (i);
901 F77_XFCN (clup1up, CLUP1UP, (m, n, F77_CMPLX_ARG (l.rwdata ()),
902 m, F77_CMPLX_ARG (r.rwdata ()), k,
903 m_ipvt.rwdata (),
904 F77_CONST_CMPLX_ARG (utmp.data ()),
905 F77_CONST_CMPLX_ARG (vtmp.data ()),
906 F77_CMPLX_ARG (w)));
907 }
908 for (F77_INT i = 0; i < m; i++) m_ipvt(i) -= 1; // decrement
909}
910
911#endif
912
913// Instantiations we need.
914
915template class lu<Matrix>;
916
917template class lu<FloatMatrix>;
918
919template class lu<ComplexMatrix>;
920
921template class lu<FloatComplexMatrix>;
922
923OCTAVE_END_NAMESPACE(math)
924OCTAVE_END_NAMESPACE(octave)
Array
N Dimensional Array with copy-on-write semantics.
Definition Array-base.h:130
Array::xelem
T & xelem(octave_idx_type n)
Size of the specified dimension.
Definition Array-base.h:547
Array::rows
octave_idx_type rows() const
Definition Array-base.h:485
Array::columns
octave_idx_type columns() const
Definition Array-base.h:497
Array::data
const T * data() const
Size of the specified dimension.
Definition Array-base.h:687
Array::rwdata
T * rwdata()
Size of the specified dimension.
Definition Array-base.cc:1790
Array::numel
octave_idx_type numel() const
Number of elements in the array.
Definition Array-base.h:440
ComplexMatrix::column
ComplexColumnVector column(octave_idx_type i) const
Definition CMatrix.cc:709
FloatComplexMatrix::column
FloatComplexColumnVector column(octave_idx_type i) const
Definition fCMatrix.cc:712
FloatMatrix::column
FloatColumnVector column(octave_idx_type i) const
Definition fMatrix.cc:429
Matrix::column
ColumnVector column(octave_idx_type i) const
Definition dMatrix.cc:423
dim_vector
Vector representing the dimensions (size) of an Array.
Definition dim-vector.h:92
lu
Definition lu.h:40
lu::packed
bool packed() const
Definition lu.cc:67
lu::lu
lu()
Definition lu.h:46
lu::update_piv
void update_piv(const VT &u, const VT &v)
lu::U
T U() const
Definition lu.cc:122
lu::unpack
void unpack()
Definition lu.cc:74
lu::P_vec
ColumnVector P_vec() const
Definition lu.cc:195
lu::regular
bool regular() const
Definition lu.cc:211
lu::update
void update(const VT &u, const VT &v)
lu::ELT_T
T::element_type ELT_T
Definition lu.h:44
lu::P
PermMatrix P() const
Definition lu.cc:188
lu::getp
Array< octave_idx_type > getp() const
Definition lu.cc:157
lu::L
T L() const
Definition lu.cc:95
lu::Y
T Y() const
Definition lu.cc:146
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.in.cc:55
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