GNU Octave 10.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-2025 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
47OCTAVE_BEGIN_NAMESPACE(octave)
48
49OCTAVE_BEGIN_NAMESPACE(math)
50
51// FIXME: PermMatrix::col_perm_vec returns Array<octave_idx_type>
52// but m_ipvt is an Array<octave_f77_int_type>. This could cause
53// trouble for large arrays if octave_f77_int_type is 32-bits but
54// octave_idx_type is 64. Since this constructor is called from
55// Fluupdate, it could be given values that are out of range. We
56// should ensure that the values are within range here.
57
58template <typename T>
59lu<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
66template <typename T>
67bool
68lu<T>::packed () const
69{
70 return m_L.dims () == dim_vector ();
71}
72
73template <typename T>
74void
75lu<T>::unpack ()
76{
77 if (packed ())
78 {
79 m_L = L ();
80 m_a_fact = U (); // FIXME: sub-optimal
81
82 // FIXME: getp returns Array<octave_idx_type> but m_ipvt is
83 // Array<octave_f77_int_type>. However, getp produces its
84 // result from a valid m_ipvt array so validation should not be
85 // necessary. OTOH, it might be better to have a version of
86 // getp that doesn't cause us to convert from
87 // Array<octave_f77_int_type> to Array<octave_idx_type> and
88 // back again.
89
90 m_ipvt = getp ();
91 }
92}
93
94template <typename T>
95T
96lu<T>::L () const
97{
98 if (packed ())
99 {
100 octave_idx_type a_nr = m_a_fact.rows ();
101 octave_idx_type a_nc = m_a_fact.columns ();
102 octave_idx_type mn = (a_nr < a_nc ? a_nr : a_nc);
103
104 T l (a_nr, mn, ELT_T (0.0));
105
106 for (octave_idx_type i = 0; i < a_nr; i++)
107 {
108 if (i < a_nc)
109 l.xelem (i, i) = 1.0;
110
111 for (octave_idx_type j = 0; j < (i < a_nc ? i : a_nc); j++)
112 l.xelem (i, j) = m_a_fact.xelem (i, j);
113 }
114
115 return l;
116 }
117 else
118 return m_L;
119}
120
121template <typename T>
122T
123lu<T>::U () const
124{
125 if (packed ())
126 {
127 octave_idx_type a_nr = m_a_fact.rows ();
128 octave_idx_type a_nc = m_a_fact.columns ();
129 octave_idx_type mn = (a_nr < a_nc ? a_nr : a_nc);
130
131 T u (mn, a_nc, ELT_T (0.0));
132
133 for (octave_idx_type i = 0; i < mn; i++)
134 {
135 for (octave_idx_type j = i; j < a_nc; j++)
136 u.xelem (i, j) = m_a_fact.xelem (i, j);
137 }
138
139 return u;
140 }
141 else
142 return m_a_fact;
143}
144
145template <typename T>
146T
147lu<T>::Y () const
148{
149 if (! packed ())
150 (*current_liboctave_error_handler)
151 ("lu: Y () not implemented for unpacked form");
152
153 return m_a_fact;
154}
155
156template <typename T>
157Array<octave_idx_type>
158lu<T>::getp () const
159{
160 if (packed ())
161 {
162 octave_idx_type a_nr = m_a_fact.rows ();
163
164 Array<octave_idx_type> pvt (dim_vector (a_nr, 1));
165
166 for (octave_idx_type i = 0; i < a_nr; i++)
167 pvt.xelem (i) = i;
168
169 for (octave_idx_type i = 0; i < m_ipvt.numel (); i++)
170 {
171 octave_idx_type k = m_ipvt.xelem (i);
172
173 if (k != i)
174 {
175 octave_idx_type tmp = pvt.xelem (k);
176 pvt.xelem (k) = pvt.xelem (i);
177 pvt.xelem (i) = tmp;
178 }
179 }
180
181 return pvt;
182 }
183 else
184 return m_ipvt;
185}
186
187template <typename T>
188PermMatrix
189lu<T>::P () const
190{
191 return PermMatrix (getp (), false);
192}
193
194template <typename T>
195ColumnVector
196lu<T>::P_vec () const
197{
198 octave_idx_type a_nr = m_a_fact.rows ();
199
200 ColumnVector p (a_nr);
201
202 Array<octave_idx_type> pvt = getp ();
203
204 for (octave_idx_type i = 0; i < a_nr; i++)
205 p.xelem (i) = static_cast<double> (pvt.xelem (i) + 1);
206
207 return p;
208}
209
210template <typename T>
211bool
212lu<T>::regular () const
213{
214 bool retval = true;
215
216 octave_idx_type k = std::min (m_a_fact.rows (), m_a_fact.columns ());
217
218 for (octave_idx_type i = 0; i < k; i++)
219 {
220 if (m_a_fact(i, i) == ELT_T ())
221 {
222 retval = false;
223 break;
224 }
225 }
226
227 return retval;
228}
229
230#if ! defined (HAVE_QRUPDATE_LUU)
231
232template <typename T>
233void
234lu<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
241template <typename T>
242void
243lu<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
250template <typename T>
251void
252lu<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
259template <typename T>
260void
261lu<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
272template <>
273OCTAVE_API
274lu<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.rwdata ();
282
283 m_a_fact = a;
284 double *tmp_data = m_a_fact.rwdata ();
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
296template <>
297OCTAVE_API void
298lu<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.rwdata (), m, r.rwdata (),
319 k, utmp.rwdata (), vtmp.rwdata ()));
320}
321
322template <>
323OCTAVE_API void
324lu<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 (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.rwdata (),
350 m, r.rwdata (), k,
351 utmp.rwdata (), vtmp.rwdata ()));
352 }
353}
354
355template <>
356OCTAVE_API void
357lu<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.rwdata (),
380 m, r.rwdata (), k,
381 m_ipvt.rwdata (),
382 utmp.data (), vtmp.data (), w));
383 for (F77_INT i = 0; i < m; i++) m_ipvt(i) -= 1; // decrement
384}
385
386template <>
387OCTAVE_API void
388lu<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 (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.rwdata (),
416 m, r.rwdata (), k,
417 m_ipvt.rwdata (),
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
425template <>
426OCTAVE_API
427lu<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.rwdata ();
435
436 m_a_fact = a;
437 float *tmp_data = m_a_fact.rwdata ();
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
449template <>
450OCTAVE_API void
451lu<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.rwdata (),
473 m, r.rwdata (), k,
474 utmp.rwdata (), vtmp.rwdata ()));
475}
476
477template <>
478OCTAVE_API void
479lu<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 (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.rwdata (),
505 m, r.rwdata (), k,
506 utmp.rwdata (), vtmp.rwdata ()));
507 }
508}
509
510template <>
511OCTAVE_API void
512lu<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.rwdata (),
536 m, r.rwdata (), k,
537 m_ipvt.rwdata (),
538 utmp.data (), vtmp.data (), w));
539 for (F77_INT i = 0; i < m; i++) m_ipvt(i) -= 1; // decrement
540}
541
542template <>
543OCTAVE_API void
544lu<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 (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.rwdata (),
572 m, r.rwdata (), k,
573 m_ipvt.rwdata (),
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
581template <>
582OCTAVE_API
583lu<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.rwdata ();
591
592 m_a_fact = a;
593 Complex *tmp_data = m_a_fact.rwdata ();
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
606template <>
607OCTAVE_API void
608lu<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.rwdata ()), m,
630 F77_DBLE_CMPLX_ARG (r.rwdata ()), k,
631 F77_DBLE_CMPLX_ARG (utmp.rwdata ()),
632 F77_DBLE_CMPLX_ARG (vtmp.rwdata ())));
633}
634
635template <>
636OCTAVE_API void
637lu<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 (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.rwdata ()),
664 m, F77_DBLE_CMPLX_ARG (r.rwdata ()),
665 k,
666 F77_DBLE_CMPLX_ARG (utmp.rwdata ()),
667 F77_DBLE_CMPLX_ARG (vtmp.rwdata ())));
668 }
669}
670
671template <>
672OCTAVE_API void
673lu<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.rwdata ()),
697 m, F77_DBLE_CMPLX_ARG (r.rwdata ()), k,
698 m_ipvt.rwdata (),
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
705template <>
706OCTAVE_API void
707lu<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 (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.rwdata ()),
737 m,
738 F77_DBLE_CMPLX_ARG (r.rwdata ()),
739 k, m_ipvt.rwdata (),
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
749template <>
750OCTAVE_API
751lu<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.rwdata ();
759
760 m_a_fact = a;
761 FloatComplex *tmp_data = m_a_fact.rwdata ();
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
774template <>
775OCTAVE_API void
776lu<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.rwdata ()), m,
798 F77_CMPLX_ARG (r.rwdata ()), k,
799 F77_CMPLX_ARG (utmp.rwdata ()),
800 F77_CMPLX_ARG (vtmp.rwdata ())));
801}
802
803template <>
804OCTAVE_API void
805lu<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 (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.rwdata ()),
832 m, F77_CMPLX_ARG (r.rwdata ()), k,
833 F77_CMPLX_ARG (utmp.rwdata ()),
834 F77_CMPLX_ARG (vtmp.rwdata ())));
835 }
836}
837
838template <>
839OCTAVE_API void
840lu<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.rwdata ()),
864 m, F77_CMPLX_ARG (r.rwdata ()), k,
865 m_ipvt.rwdata (),
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
872template <>
873OCTAVE_API void
874lu<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 (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.rwdata ()),
903 m, F77_CMPLX_ARG (r.rwdata ()), k,
904 m_ipvt.rwdata (),
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
916template class lu<Matrix>;
917
918template class lu<FloatMatrix>;
919
920template class lu<ComplexMatrix>;
921
922template class lu<FloatComplexMatrix>;
923
924OCTAVE_END_NAMESPACE(math)
925OCTAVE_END_NAMESPACE(octave)
Array
N Dimensional Array with copy-on-write semantics.
Definition Array.h:130
Array::xelem
T & xelem(octave_idx_type n)
Size of the specified dimension.
Definition Array.h:525
Array::rows
octave_idx_type rows() const
Definition Array.h:463
Array::columns
octave_idx_type columns() const
Definition Array.h:475
Array::data
const T * data() const
Size of the specified dimension.
Definition Array.h:665
Array::rwdata
T * rwdata()
Size of the specified dimension.
Definition Array-base.cc:1787
Array::numel
octave_idx_type numel() const
Number of elements in the array.
Definition Array.h:418
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:90
lu
Definition lu.h:41
lu::packed
bool packed() const
Definition lu.cc:68
lu::lu
lu()
Definition lu.h:47
lu::update_piv
void update_piv(const VT &u, const VT &v)
lu::U
T U() const
Definition lu.cc:123
lu::unpack
void unpack()
Definition lu.cc:75
lu::P_vec
ColumnVector P_vec() const
Definition lu.cc:196
lu::regular
bool regular() const
Definition lu.cc:212
lu::update
void update(const VT &u, const VT &v)
lu::ELT_T
T::element_type ELT_T
Definition lu.h:45
lu::P
PermMatrix P() const
Definition lu.cc:189
lu::getp
Array< octave_idx_type > getp() const
Definition lu.cc:158
lu::L
T L() const
Definition lu.cc:96
lu::Y
T Y() const
Definition lu.cc:147
octave_idx_type
OCTAVE_BEGIN_NAMESPACE
OCTAVE_BEGIN_NAMESPACE(octave) static octave_value daspk_fcn
F77_CONST_CMPLX_ARG
#define F77_CONST_CMPLX_ARG(x)
Definition f77-fcn.h:313
F77_DBLE_CMPLX_ARG
#define F77_DBLE_CMPLX_ARG(x)
Definition f77-fcn.h:316
F77_CMPLX_ARG
#define F77_CMPLX_ARG(x)
Definition f77-fcn.h:310
F77_XFCN
#define F77_XFCN(f, F, args)
Definition f77-fcn.h:45
F77_INT
octave_f77_int_type F77_INT
Definition f77-fcn.h:306
F77_CONST_DBLE_CMPLX_ARG
#define F77_CONST_DBLE_CMPLX_ARG(x)
Definition f77-fcn.h:319
OCTAVE_API
#define OCTAVE_API
Definition main.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