Mercurial > octave-nkf

comparison src/DLD-FUNCTIONS/luinc.cc @ 10154:40dfc0c99116

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
DLD-FUNCTIONS/*.cc: untabify
author John W. Eaton <jwe@octave.org>
date Wed, 20 Jan 2010 17:33:41 -0500
parents 7c02ec148a3c
children d0ce5e973937
comparison
equal deleted inserted replaced
10153:2c28f9d0360f 10154:40dfc0c99116
111 Matrix thresh; 111 Matrix thresh;
112 double droptol = -1.; 112 double droptol = -1.;
113 bool vecout; 113 bool vecout;
114 114
115 if (args(1).is_string ()) 115 if (args(1).is_string ())
116 { 116 {
117 if (args(1).string_value () == "0") 117 if (args(1).string_value () == "0")
118 zero_level = true; 118 zero_level = true;
119 else 119 else
120 error ("luinc: unrecognized string argument"); 120 error ("luinc: unrecognized string argument");
121 } 121 }
122 else if (args(1).is_map ()) 122 else if (args(1).is_map ())
123 { 123 {
124 Octave_map map = args(1).map_value (); 124 Octave_map map = args(1).map_value ();
125 125
126 if (map.contains ("droptol")) 126 if (map.contains ("droptol"))
127 droptol = map.contents ("droptol")(0).double_value (); 127 droptol = map.contents ("droptol")(0).double_value ();
128 128
129 if (map.contains ("milu")) 129 if (map.contains ("milu"))
130 { 130 {
131 double tmp = map.contents ("milu")(0).double_value (); 131 double tmp = map.contents ("milu")(0).double_value ();
132 132
133 milu = (tmp == 0. ? false : true); 133 milu = (tmp == 0. ? false : true);
134 } 134 }
135 135
136 if (map.contains ("udiag")) 136 if (map.contains ("udiag"))
137 { 137 {
138 double tmp = map.contents ("udiag")(0).double_value (); 138 double tmp = map.contents ("udiag")(0).double_value ();
139 139
140 udiag = (tmp == 0. ? false : true); 140 udiag = (tmp == 0. ? false : true);
141 } 141 }
142 142
143 if (map.contains ("thresh")) 143 if (map.contains ("thresh"))
144 { 144 {
145 thresh = map.contents ("thresh")(0).matrix_value (); 145 thresh = map.contents ("thresh")(0).matrix_value ();
146 146
147 if (thresh.nelem () == 1) 147 if (thresh.nelem () == 1)
148 { 148 {
149 thresh.resize(1,2); 149 thresh.resize(1,2);
150 thresh(1) = thresh(0); 150 thresh(1) = thresh(0);
151 } 151 }
152 else if (thresh.nelem () != 2) 152 else if (thresh.nelem () != 2)
153 error ("chol: expecting 2 element vector for thresh"); 153 error ("chol: expecting 2 element vector for thresh");
154 } 154 }
155 } 155 }
156 else 156 else
157 droptol = args(1).double_value (); 157 droptol = args(1).double_value ();
158 158
159 if (nargin == 3) 159 if (nargin == 3)
160 { 160 {
161 std::string tmp = args(2).string_value (); 161 std::string tmp = args(2).string_value ();
162 162
163 if (! error_state ) 163 if (! error_state )
164 { 164 {
165 if (tmp.compare ("vector") == 0) 165 if (tmp.compare ("vector") == 0)
166 vecout = true; 166 vecout = true;
167 else 167 else
168 error ("luinc: unrecognized string argument"); 168 error ("luinc: unrecognized string argument");
169 } 169 }
170 } 170 }
171 171
172 // FIXME Add code for zero-level factorization 172 // FIXME Add code for zero-level factorization
173 if (zero_level) 173 if (zero_level)
174 error ("luinc: zero-level factorization not implemented"); 174 error ("luinc: zero-level factorization not implemented");
175 175
176 if (!error_state) 176 if (!error_state)
177 { 177 {
178 if (args(0).type_name () == "sparse matrix") 178 if (args(0).type_name () == "sparse matrix")
179 { 179 {
180 SparseMatrix sm = args(0).sparse_matrix_value (); 180 SparseMatrix sm = args(0).sparse_matrix_value ();
181 octave_idx_type sm_nr = sm.rows (); 181 octave_idx_type sm_nr = sm.rows ();
182 octave_idx_type sm_nc = sm.cols (); 182 octave_idx_type sm_nc = sm.cols ();
183 ColumnVector Qinit (sm_nc); 183 ColumnVector Qinit (sm_nc);
184 184
185 for (octave_idx_type i = 0; i < sm_nc; i++) 185 for (octave_idx_type i = 0; i < sm_nc; i++)
186 Qinit (i) = i; 186 Qinit (i) = i;
187 187
188 if (! error_state) 188 if (! error_state)
189 { 189 {
190 switch (nargout) 190 switch (nargout)
191 { 191 {
192 case 0: 192 case 0:
193 case 1: 193 case 1:
194 case 2: 194 case 2:
195 { 195 {
196 SparseLU fact (sm, Qinit, thresh, false, true, droptol, 196 SparseLU fact (sm, Qinit, thresh, false, true, droptol,
197 milu, udiag); 197 milu, udiag);
198 198
199 if (! error_state) 199 if (! error_state)
200 { 200 {
201 SparseMatrix P = fact.Pr (); 201 SparseMatrix P = fact.Pr ();
202 SparseMatrix L = P.transpose () * fact.L (); 202 SparseMatrix L = P.transpose () * fact.L ();
203 retval(1) = octave_value (fact.U (), 203 retval(1) = octave_value (fact.U (),
204 MatrixType (MatrixType::Upper)); 204 MatrixType (MatrixType::Upper));
205 retval(0) = octave_value (L, MatrixType 205 retval(0) = octave_value (L, MatrixType
206 (MatrixType::Permuted_Lower, 206 (MatrixType::Permuted_Lower,
207 sm_nr, fact.row_perm ())); 207 sm_nr, fact.row_perm ()));
208 } 208 }
209 } 209 }
210 break; 210 break;
211 211
212 case 3: 212 case 3:
213 { 213 {
214 SparseLU fact (sm, Qinit, thresh, false, true, droptol, 214 SparseLU fact (sm, Qinit, thresh, false, true, droptol,
215 milu, udiag); 215 milu, udiag);
216 216
217 if (! error_state) 217 if (! error_state)
218 { 218 {
219 if (vecout) 219 if (vecout)
220 retval(2) = fact.Pr_vec (); 220 retval(2) = fact.Pr_vec ();
221 else 221 else
222 retval(2) = fact.Pr (); 222 retval(2) = fact.Pr ();
223 retval(1) = octave_value (fact.U (), 223 retval(1) = octave_value (fact.U (),
224 MatrixType (MatrixType::Upper)); 224 MatrixType (MatrixType::Upper));
225 retval(0) = octave_value (fact.L (), 225 retval(0) = octave_value (fact.L (),
226 MatrixType (MatrixType::Lower)); 226 MatrixType (MatrixType::Lower));
227 } 227 }
228 } 228 }
229 break; 229 break;
230 230
231 case 4: 231 case 4:
232 default: 232 default:
233 { 233 {
234 SparseLU fact (sm, Qinit, thresh, false, false, droptol, 234 SparseLU fact (sm, Qinit, thresh, false, false, droptol,
235 milu, udiag); 235 milu, udiag);
236 236
237 if (! error_state) 237 if (! error_state)
238 { 238 {
239 if (vecout) 239 if (vecout)
240 { 240 {
241 retval(3) = fact.Pc_vec (); 241 retval(3) = fact.Pc_vec ();
242 retval(2) = fact.Pr_vec (); 242 retval(2) = fact.Pr_vec ();
243 } 243 }
244 else 244 else
245 { 245 {
246 retval(3) = fact.Pc (); 246 retval(3) = fact.Pc ();
247 retval(2) = fact.Pr (); 247 retval(2) = fact.Pr ();
248 } 248 }
249 retval(1) = octave_value (fact.U (), 249 retval(1) = octave_value (fact.U (),
250 MatrixType (MatrixType::Upper)); 250 MatrixType (MatrixType::Upper));
251 retval(0) = octave_value (fact.L (), 251 retval(0) = octave_value (fact.L (),
252 MatrixType (MatrixType::Lower)); 252 MatrixType (MatrixType::Lower));
253 } 253 }
254 } 254 }
255 break; 255 break;
256 } 256 }
257 } 257 }
258 } 258 }
259 else if (args(0).type_name () == "sparse complex matrix") 259 else if (args(0).type_name () == "sparse complex matrix")
260 { 260 {
261 SparseComplexMatrix sm = 261 SparseComplexMatrix sm =
262 args(0).sparse_complex_matrix_value (); 262 args(0).sparse_complex_matrix_value ();
263 octave_idx_type sm_nr = sm.rows (); 263 octave_idx_type sm_nr = sm.rows ();
264 octave_idx_type sm_nc = sm.cols (); 264 octave_idx_type sm_nc = sm.cols ();
265 ColumnVector Qinit (sm_nc); 265 ColumnVector Qinit (sm_nc);
266 266
267 for (octave_idx_type i = 0; i < sm_nc; i++) 267 for (octave_idx_type i = 0; i < sm_nc; i++)
268 Qinit (i) = i; 268 Qinit (i) = i;
269 269
270 if (! error_state) 270 if (! error_state)
271 { 271 {
272 switch (nargout) 272 switch (nargout)
273 { 273 {
274 case 0: 274 case 0:
275 case 1: 275 case 1:
276 case 2: 276 case 2:
277 { 277 {
278 SparseComplexLU fact (sm, Qinit, thresh, false, true, 278 SparseComplexLU fact (sm, Qinit, thresh, false, true,
279 droptol, milu, udiag); 279 droptol, milu, udiag);
280 280
281 281
282 if (! error_state) 282 if (! error_state)
283 { 283 {
284 SparseMatrix P = fact.Pr (); 284 SparseMatrix P = fact.Pr ();
285 SparseComplexMatrix L = P.transpose () * fact.L (); 285 SparseComplexMatrix L = P.transpose () * fact.L ();
286 retval(1) = octave_value (fact.U (), 286 retval(1) = octave_value (fact.U (),
287 MatrixType (MatrixType::Upper)); 287 MatrixType (MatrixType::Upper));
288 retval(0) = octave_value (L, MatrixType 288 retval(0) = octave_value (L, MatrixType
289 (MatrixType::Permuted_Lower, 289 (MatrixType::Permuted_Lower,
290 sm_nr, fact.row_perm ())); 290 sm_nr, fact.row_perm ()));
291 } 291 }
292 } 292 }
293 break; 293 break;
294 294
295 case 3: 295 case 3:
296 { 296 {
297 SparseComplexLU fact (sm, Qinit, thresh, false, true, 297 SparseComplexLU fact (sm, Qinit, thresh, false, true,
298 droptol, milu, udiag); 298 droptol, milu, udiag);
299 299
300 if (! error_state) 300 if (! error_state)
301 { 301 {
302 if (vecout) 302 if (vecout)
303 retval(2) = fact.Pr_vec (); 303 retval(2) = fact.Pr_vec ();
304 else 304 else
305 retval(2) = fact.Pr (); 305 retval(2) = fact.Pr ();
306 retval(1) = octave_value (fact.U (), 306 retval(1) = octave_value (fact.U (),
307 MatrixType (MatrixType::Upper)); 307 MatrixType (MatrixType::Upper));
308 retval(0) = octave_value (fact.L (), 308 retval(0) = octave_value (fact.L (),
309 MatrixType (MatrixType::Lower)); 309 MatrixType (MatrixType::Lower));
310 } 310 }
311 } 311 }
312 break; 312 break;
313 313
314 case 4: 314 case 4:
315 default: 315 default:
316 { 316 {
317 SparseComplexLU fact (sm, Qinit, thresh, false, false, 317 SparseComplexLU fact (sm, Qinit, thresh, false, false,
318 droptol, milu, udiag); 318 droptol, milu, udiag);
319 319
320 if (! error_state) 320 if (! error_state)
321 { 321 {
322 if (vecout) 322 if (vecout)
323 { 323 {
324 retval(3) = fact.Pc_vec (); 324 retval(3) = fact.Pc_vec ();
325 retval(2) = fact.Pr_vec (); 325 retval(2) = fact.Pr_vec ();
326 } 326 }
327 else 327 else
328 { 328 {
329 retval(3) = fact.Pc (); 329 retval(3) = fact.Pc ();
330 retval(2) = fact.Pr (); 330 retval(2) = fact.Pr ();
331 } 331 }
332 retval(1) = octave_value (fact.U (), 332 retval(1) = octave_value (fact.U (),
333 MatrixType (MatrixType::Upper)); 333 MatrixType (MatrixType::Upper));
334 retval(0) = octave_value (fact.L (), 334 retval(0) = octave_value (fact.L (),
335 MatrixType (MatrixType::Lower)); 335 MatrixType (MatrixType::Lower));
336 } 336 }
337 } 337 }
338 break; 338 break;
339 } 339 }
340 } 340 }
341 } 341 }
342 else 342 else
343 error ("luinc: first argument must be sparse"); 343 error ("luinc: first argument must be sparse");
344 } 344 }
345 } 345 }
346 346
347 return retval; 347 return retval;
348 } 348 }
349 349