Mercurial > octave-nkf

comparison src/DLD-FUNCTIONS/luinc.cc @ 5322:22994a5730f9

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 2005-04-29 13:04:24 by dbateman]
author dbateman
date Fri, 29 Apr 2005 13:04:25 +0000
parents 4c8a2e4e0717
children 9761b7d24e9e
comparison
equal deleted inserted replaced
5321:84b72a402b86 5322:22994a5730f9
28 #include "gripes.h" 28 #include "gripes.h"
29 #include "oct-obj.h" 29 #include "oct-obj.h"
30 #include "utils.h" 30 #include "utils.h"
31 #include "oct-map.h" 31 #include "oct-map.h"
32 32
33 #include "SparseType.h"
33 #include "SparseCmplxLU.h" 34 #include "SparseCmplxLU.h"
34 #include "SparsedbleLU.h" 35 #include "SparsedbleLU.h"
35 #include "ov-re-sparse.h" 36 #include "ov-re-sparse.h"
36 #include "ov-cx-sparse.h" 37 #include "ov-cx-sparse.h"
37 38
144 if (zero_level) 145 if (zero_level)
145 error ("luinc: zero-level factorization not implemented"); 146 error ("luinc: zero-level factorization not implemented");
146 147
147 if (!error_state) 148 if (!error_state)
148 { 149 {
149 if (args(0).class_name () == "sparse") 150 if (args(0).type_name () == "sparse matrix")
150 { 151 {
151 SparseMatrix sm = args(0).sparse_matrix_value (); 152 SparseMatrix sm = args(0).sparse_matrix_value ();
153 octave_idx_type sm_nr = sm.rows ();
152 octave_idx_type sm_nc = sm.cols (); 154 octave_idx_type sm_nc = sm.cols ();
153 ColumnVector Qinit (sm_nc); 155 ColumnVector Qinit (sm_nc);
154 156
155 for (octave_idx_type i = 0; i < sm_nc; i++) 157 for (octave_idx_type i = 0; i < sm_nc; i++)
156 Qinit (i) = i; 158 Qinit (i) = i;
166 SparseLU fact (sm, Qinit, thresh, true, droptol, 168 SparseLU fact (sm, Qinit, thresh, true, droptol,
167 milu, udiag); 169 milu, udiag);
168 170
169 SparseMatrix P = fact.Pr (); 171 SparseMatrix P = fact.Pr ();
170 SparseMatrix L = P.transpose () * fact.L (); 172 SparseMatrix L = P.transpose () * fact.L ();
171 retval(1) = fact.U (); 173 retval(1) = octave_value (fact.U (),
172 retval(0) = L; 174 SparseType (SparseType::Upper));
175 retval(0) = octave_value (L, SparseType
176 (SparseType::Permuted_Lower,
177 sm_nr, fact.row_perm ()));
173 } 178 }
174 break; 179 break;
175 180
176 case 3: 181 case 3:
177 { 182 {
178 SparseLU fact (sm, Qinit, thresh, true, droptol, 183 SparseLU fact (sm, Qinit, thresh, true, droptol,
179 milu, udiag); 184 milu, udiag);
180 185
181 retval(2) = fact.Pr (); 186 retval(2) = fact.Pr ();
182 retval(1) = fact.U (); 187 retval(1) = octave_value (fact.U (),
183 retval(0) = fact.L (); 188 SparseType (SparseType::Upper));
189 retval(0) = octave_value (fact.L (),
190 SparseType (SparseType::Lower));
184 } 191 }
185 break; 192 break;
186 193
187 case 4: 194 case 4:
188 default: 195 default:
190 SparseLU fact (sm, Qinit, thresh, false, droptol, 197 SparseLU fact (sm, Qinit, thresh, false, droptol,
191 milu, udiag); 198 milu, udiag);
192 199
193 retval(3) = fact.Pc (); 200 retval(3) = fact.Pc ();
194 retval(2) = fact.Pr (); 201 retval(2) = fact.Pr ();
195 retval(1) = fact.U (); 202 retval(1) = octave_value (fact.U (),
196 retval(0) = fact.L (); 203 SparseType (SparseType::Upper));
204 retval(0) = octave_value (fact.L (),
205 SparseType (SparseType::Lower));
197 } 206 }
198 break; 207 break;
199 } 208 }
200 } 209 }
201 } 210 }
202 else if (args(0).class_name () == "sparse complex") 211 else if (args(0).type_name () == "sparse complex matrix")
203 { 212 {
204 SparseComplexMatrix sm = 213 SparseComplexMatrix sm =
205 args(0).sparse_complex_matrix_value (); 214 args(0).sparse_complex_matrix_value ();
215 octave_idx_type sm_nr = sm.rows ();
206 octave_idx_type sm_nc = sm.cols (); 216 octave_idx_type sm_nc = sm.cols ();
207 ColumnVector Qinit (sm_nc); 217 ColumnVector Qinit (sm_nc);
208 218
209 for (octave_idx_type i = 0; i < sm_nc; i++) 219 for (octave_idx_type i = 0; i < sm_nc; i++)
210 Qinit (i) = i; 220 Qinit (i) = i;
220 SparseComplexLU fact (sm, Qinit, thresh, true, 230 SparseComplexLU fact (sm, Qinit, thresh, true,
221 droptol, milu, udiag); 231 droptol, milu, udiag);
222 232
223 SparseMatrix P = fact.Pr (); 233 SparseMatrix P = fact.Pr ();
224 SparseComplexMatrix L = P.transpose () * fact.L (); 234 SparseComplexMatrix L = P.transpose () * fact.L ();
225 retval(1) = fact.U (); 235 retval(1) = octave_value (fact.U (),
226 retval(0) = L; 236 SparseType (SparseType::Upper));
237 retval(0) = octave_value (L, SparseType
238 (SparseType::Permuted_Lower,
239 sm_nr, fact.row_perm ()));
227 } 240 }
228 break; 241 break;
229 242
230 case 3: 243 case 3:
231 { 244 {
232 SparseComplexLU fact (sm, Qinit, thresh, true, 245 SparseComplexLU fact (sm, Qinit, thresh, true,
233 droptol, milu, udiag); 246 droptol, milu, udiag);
234 247
235 retval(2) = fact.Pr (); 248 retval(2) = fact.Pr ();
236 retval(1) = fact.U (); 249 retval(1) = octave_value (fact.U (),
237 retval(0) = fact.L (); 250 SparseType (SparseType::Upper));
251 retval(0) = octave_value (fact.L (),
252 SparseType (SparseType::Lower));
238 } 253 }
239 break; 254 break;
240 255
241 case 4: 256 case 4:
242 default: 257 default:
244 SparseComplexLU fact (sm, Qinit, thresh, false, 259 SparseComplexLU fact (sm, Qinit, thresh, false,
245 droptol, milu, udiag); 260 droptol, milu, udiag);
246 261
247 retval(3) = fact.Pc (); 262 retval(3) = fact.Pc ();
248 retval(2) = fact.Pr (); 263 retval(2) = fact.Pr ();
249 retval(1) = fact.U (); 264 retval(1) = octave_value (fact.U (),
250 retval(0) = fact.L (); 265 SparseType (SparseType::Upper));
266 retval(0) = octave_value (fact.L (),
267 SparseType (SparseType::Lower));
251 } 268 }
252 break; 269 break;
253 } 270 }
254 } 271 }
255 } 272 }