GNU Octave  3.8.0
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
CmplxLU.cc
Go to the documentation of this file.
1 /*
2 
3 Copyright (C) 1994-2013 John W. Eaton
4 Copyright (C) 2009 VZLU Prague, a.s.
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 #ifdef HAVE_CONFIG_H
25 #include <config.h>
26 #endif
27 
28 #include "CmplxLU.h"
29 #include "f77-fcn.h"
30 #include "lo-error.h"
31 #include "oct-locbuf.h"
32 #include "CColVector.h"
33 
34 // Instantiate the base LU class for the types we need.
35 
36 #include "base-lu.h"
37 #include "base-lu.cc"
38 
39 template class base_lu <ComplexMatrix>;
40 
41 // Define the constructor for this particular derivation.
42 
43 extern "C"
44 {
45  F77_RET_T
46  F77_FUNC (zgetrf, ZGETRF) (const octave_idx_type&, const octave_idx_type&,
47  Complex*, const octave_idx_type&,
48  octave_idx_type*, octave_idx_type&);
49 
50 #ifdef HAVE_QRUPDATE_LUU
51  F77_RET_T
52  F77_FUNC (zlu1up, ZLU1UP) (const octave_idx_type&, const octave_idx_type&,
53  Complex *, const octave_idx_type&,
54  Complex *, const octave_idx_type&,
55  Complex *, Complex *);
56 
57  F77_RET_T
58  F77_FUNC (zlup1up, ZLUP1UP) (const octave_idx_type&, const octave_idx_type&,
59  Complex *, const octave_idx_type&,
60  Complex *, const octave_idx_type&,
61  octave_idx_type *, const Complex *,
62  const Complex *, Complex *);
63 #endif
64 }
65 
67 {
68  octave_idx_type a_nr = a.rows ();
69  octave_idx_type a_nc = a.cols ();
70  octave_idx_type mn = (a_nr < a_nc ? a_nr : a_nc);
71 
72  ipvt.resize (dim_vector (mn, 1));
73  octave_idx_type *pipvt = ipvt.fortran_vec ();
74 
75  a_fact = a;
76  Complex *tmp_data = a_fact.fortran_vec ();
77 
78  octave_idx_type info = 0;
79 
80  F77_XFCN (zgetrf, ZGETRF, (a_nr, a_nc, tmp_data, a_nr, pipvt, info));
81 
82  for (octave_idx_type i = 0; i < mn; i++)
83  pipvt[i] -= 1;
84 }
85 
86 #ifdef HAVE_QRUPDATE_LUU
87 
89  const ComplexColumnVector& v)
90 {
91  if (packed ())
92  unpack ();
93 
94  ComplexMatrix& l = l_fact;
95  ComplexMatrix& r = a_fact;
96 
97  octave_idx_type m = l.rows ();
98  octave_idx_type n = r.columns ();
99  octave_idx_type k = l.columns ();
100 
101  if (u.length () == m && v.length () == n)
102  {
103  ComplexColumnVector utmp = u, vtmp = v;
104  F77_XFCN (zlu1up, ZLU1UP, (m, n, l.fortran_vec (), m, r.fortran_vec (), k,
105  utmp.fortran_vec (), vtmp.fortran_vec ()));
106  }
107  else
108  (*current_liboctave_error_handler) ("luupdate: dimensions mismatch");
109 }
110 
112 {
113  if (packed ())
114  unpack ();
115 
116  ComplexMatrix& l = l_fact;
117  ComplexMatrix& r = a_fact;
118 
119  octave_idx_type m = l.rows ();
120  octave_idx_type n = r.columns ();
121  octave_idx_type k = l.columns ();
122 
123  if (u.rows () == m && v.rows () == n && u.cols () == v.cols ())
124  {
125  for (volatile octave_idx_type i = 0; i < u.cols (); i++)
126  {
127  ComplexColumnVector utmp = u.column (i), vtmp = v.column (i);
128  F77_XFCN (zlu1up, ZLU1UP, (m, n, l.fortran_vec (),
129  m, r.fortran_vec (), k,
130  utmp.fortran_vec (), vtmp.fortran_vec ()));
131  }
132  }
133  else
134  (*current_liboctave_error_handler) ("luupdate: dimensions mismatch");
135 }
136 
138  const ComplexColumnVector& v)
139 {
140  if (packed ())
141  unpack ();
142 
143  ComplexMatrix& l = l_fact;
144  ComplexMatrix& r = a_fact;
145 
146  octave_idx_type m = l.rows ();
147  octave_idx_type n = r.columns ();
148  octave_idx_type k = l.columns ();
149 
150  if (u.length () == m && v.length () == n)
151  {
152  ComplexColumnVector utmp = u, vtmp = v;
154  for (octave_idx_type i = 0; i < m; i++) ipvt(i) += 1; // increment
155  F77_XFCN (zlup1up, ZLUP1UP, (m, n, l.fortran_vec (),
156  m, r.fortran_vec (), k,
157  ipvt.fortran_vec (),
158  utmp.data (), vtmp.data (), w));
159  for (octave_idx_type i = 0; i < m; i++) ipvt(i) -= 1; // decrement
160  }
161  else
162  (*current_liboctave_error_handler) ("luupdate: dimensions mismatch");
163 }
164 
166 {
167  if (packed ())
168  unpack ();
169 
170  ComplexMatrix& l = l_fact;
171  ComplexMatrix& r = a_fact;
172 
173  octave_idx_type m = l.rows ();
174  octave_idx_type n = r.columns ();
175  octave_idx_type k = l.columns ();
176 
177  if (u.rows () == m && v.rows () == n && u.cols () == v.cols ())
178  {
180  for (octave_idx_type i = 0; i < m; i++) ipvt(i) += 1; // increment
181  for (volatile octave_idx_type i = 0; i < u.cols (); i++)
182  {
183  ComplexColumnVector utmp = u.column (i), vtmp = v.column (i);
184  F77_XFCN (zlup1up, ZLUP1UP, (m, n, l.fortran_vec (),
185  m, r.fortran_vec (), k,
186  ipvt.fortran_vec (),
187  utmp.data (), vtmp.data (), w));
188  }
189  for (octave_idx_type i = 0; i < m; i++) ipvt(i) -= 1; // decrement
190  }
191  else
192  (*current_liboctave_error_handler) ("luupdate: dimensions mismatch");
193 }
194 
195 #else
196 
198 {
199  (*current_liboctave_error_handler)
200  ("luupdate: not available in this version");
201 }
202 
203 void ComplexLU::update (const ComplexMatrix&, const ComplexMatrix&)
204 {
205  (*current_liboctave_error_handler)
206  ("luupdate: not available in this version");
207 }
208 
210  const ComplexColumnVector&)
211 {
212  (*current_liboctave_error_handler)
213  ("luupdate: not available in this version");
214 }
215 
217 {
218  (*current_liboctave_error_handler)
219  ("luupdate: not available in this version");
220 }
221 
222 #endif