Mercurial > jwe > octave

comparison liboctave/numeric/fCmplxCHOL.cc @ 21136:7cac4e7458f2

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
maint: clean up code around calls to current_liboctave_error_handler. Remove statements after call to handler that are no longer reachable. Place input validation first and immediately call handler if necessary. Change if/error_handler/else to if/error_handler and re-indent code. * Array-util.cc, Array.cc, CColVector.cc, CDiagMatrix.cc, CMatrix.cc, CNDArray.cc, CRowVector.cc, CSparse.cc, DiagArray2.cc, MArray.cc, PermMatrix.cc, Sparse.cc, Sparse.h, chMatrix.cc, chNDArray.cc, dColVector.cc, dDiagMatrix.cc, dMatrix.cc, dNDArray.cc, dRowVector.cc, dSparse.cc, fCColVector.cc, fCDiagMatrix.cc, fCMatrix.cc, fCNDArray.cc, fCRowVector.cc, fColVector.cc, fDiagMatrix.cc, fMatrix.cc, fNDArray.cc, fRowVector.cc, idx-vector.cc, CmplxAEPBAL.cc, CmplxCHOL.cc, CmplxGEPBAL.cc, CmplxHESS.cc, CmplxLU.cc, CmplxQR.cc, CmplxSCHUR.cc, CmplxSVD.cc, DASPK.cc, EIG.cc, LSODE.cc, Quad.cc, SparseCmplxCHOL.cc, SparseCmplxLU.cc, SparseCmplxQR.cc, SparseQR.cc, SparsedbleCHOL.cc, SparsedbleLU.cc, base-lu.cc, bsxfun-defs.cc, dbleAEPBAL.cc, dbleCHOL.cc, dbleGEPBAL.cc, dbleHESS.cc, dbleLU.cc, dbleQR.cc, dbleSCHUR.cc, dbleSVD.cc, eigs-base.cc, fCmplxAEPBAL.cc, fCmplxCHOL.cc, fCmplxLU.cc, fCmplxQR.cc, fCmplxSCHUR.cc, fEIG.cc, floatAEPBAL.cc, floatCHOL.cc, floatGEPBAL.cc, floatHESS.cc, floatLU.cc, floatQR.cc, floatSCHUR.cc, floatSVD.cc, lo-specfun.cc, oct-fftw.cc, oct-rand.cc, oct-spparms.cc, sparse-base-chol.cc, sparse-dmsolve.cc, file-ops.cc, lo-sysdep.cc, mach-info.cc, oct-env.cc, oct-syscalls.cc, cmd-edit.cc, cmd-hist.cc, data-conv.cc, lo-ieee.cc, lo-regexp.cc, oct-base64.cc, oct-shlib.cc, pathsearch.cc, singleton-cleanup.cc, sparse-util.cc, unwind-prot.cc: Remove statements after call to handler that are no longer reachable. Place input validation first and immediately call handler if necessary. Change if/error_handler/else to if/error_handler and re-indent code.
author Rik <rik@octave.org>
date Sat, 23 Jan 2016 13:52:03 -0800
parents ab705b42cfd8
children f7121e111991
comparison
equal deleted inserted replaced
21135:95da3bc8a281 21136:7cac4e7458f2
90 { 90 {
91 octave_idx_type a_nr = a.rows (); 91 octave_idx_type a_nr = a.rows ();
92 octave_idx_type a_nc = a.cols (); 92 octave_idx_type a_nc = a.cols ();
93 93
94 if (a_nr != a_nc) 94 if (a_nr != a_nc)
95 { 95 (*current_liboctave_error_handler)
96 (*current_liboctave_error_handler) 96 ("FloatComplexCHOL requires square matrix");
97 ("FloatComplexCHOL requires square matrix");
98 return -1;
99 }
100 97
101 octave_idx_type n = a_nc; 98 octave_idx_type n = a_nc;
102 octave_idx_type info; 99 octave_idx_type info;
103 100
104 is_upper = upper; 101 is_upper = upper;
163 FloatComplexMatrix retval; 160 FloatComplexMatrix retval;
164 161
165 octave_idx_type r_nr = r.rows (); 162 octave_idx_type r_nr = r.rows ();
166 octave_idx_type r_nc = r.cols (); 163 octave_idx_type r_nc = r.cols ();
167 164
168 if (r_nr == r_nc) 165 if (r_nr != r_nc)
169 { 166 (*current_liboctave_error_handler) ("chol2inv requires square matrix");
170 octave_idx_type n = r_nc; 167
171 octave_idx_type info; 168 octave_idx_type n = r_nc;
172 169 octave_idx_type info;
173 FloatComplexMatrix tmp = r; 170
174 171 FloatComplexMatrix tmp = r;
172
173 if (is_upper)
174 F77_XFCN (cpotri, CPOTRI, (F77_CONST_CHAR_ARG2 ("U", 1), n,
175 tmp.fortran_vec (), n, info
176 F77_CHAR_ARG_LEN (1)));
177 else
178 F77_XFCN (cpotri, CPOTRI, (F77_CONST_CHAR_ARG2 ("L", 1), n,
179 tmp.fortran_vec (), n, info
180 F77_CHAR_ARG_LEN (1)));
181
182 // If someone thinks of a more graceful way of doing this (or
183 // faster for that matter :-)), please let me know!
184
185 if (n > 1)
186 {
175 if (is_upper) 187 if (is_upper)
176 F77_XFCN (cpotri, CPOTRI, (F77_CONST_CHAR_ARG2 ("U", 1), n, 188 for (octave_idx_type j = 0; j < r_nc; j++)
177 tmp.fortran_vec (), n, info 189 for (octave_idx_type i = j+1; i < r_nr; i++)
178 F77_CHAR_ARG_LEN (1))); 190 tmp.xelem (i, j) = tmp.xelem (j, i);
179 else 191 else
180 F77_XFCN (cpotri, CPOTRI, (F77_CONST_CHAR_ARG2 ("L", 1), n, 192 for (octave_idx_type j = 0; j < r_nc; j++)
181 tmp.fortran_vec (), n, info 193 for (octave_idx_type i = j+1; i < r_nr; i++)
182 F77_CHAR_ARG_LEN (1))); 194 tmp.xelem (j, i) = tmp.xelem (i, j);
183 195 }
184 // If someone thinks of a more graceful way of doing this (or 196
185 // faster for that matter :-)), please let me know! 197 retval = tmp;
186
187 if (n > 1)
188 {
189 if (is_upper)
190 for (octave_idx_type j = 0; j < r_nc; j++)
191 for (octave_idx_type i = j+1; i < r_nr; i++)
192 tmp.xelem (i, j) = tmp.xelem (j, i);
193 else
194 for (octave_idx_type j = 0; j < r_nc; j++)
195 for (octave_idx_type i = j+1; i < r_nr; i++)
196 tmp.xelem (j, i) = tmp.xelem (i, j);
197 }
198
199 retval = tmp;
200 }
201 else
202 (*current_liboctave_error_handler) ("chol2inv requires square matrix");
203 198
204 return retval; 199 return retval;
205 } 200 }
206 201
207 // Compute the inverse of a matrix using the Cholesky factorization. 202 // Compute the inverse of a matrix using the Cholesky factorization.
212 } 207 }
213 208
214 void 209 void
215 FloatComplexCHOL::set (const FloatComplexMatrix& R) 210 FloatComplexCHOL::set (const FloatComplexMatrix& R)
216 { 211 {
217 if (R.is_square ()) 212 if (! R.is_square ())
218 chol_mat = R;
219 else
220 (*current_liboctave_error_handler) ("CHOL requires square matrix"); 213 (*current_liboctave_error_handler) ("CHOL requires square matrix");
214
215 chol_mat = R;
221 } 216 }
222 217
223 #ifdef HAVE_QRUPDATE 218 #ifdef HAVE_QRUPDATE
224 219
225 void 220 void
226 FloatComplexCHOL::update (const FloatComplexColumnVector& u) 221 FloatComplexCHOL::update (const FloatComplexColumnVector& u)
227 { 222 {
228 octave_idx_type n = chol_mat.rows (); 223 octave_idx_type n = chol_mat.rows ();
229 224
230 if (u.numel () == n) 225 if (u.numel () != n)
231 {
232 FloatComplexColumnVector utmp = u;
233
234 OCTAVE_LOCAL_BUFFER (float, rw, n);
235
236 F77_XFCN (cch1up, CCH1UP, (n, chol_mat.fortran_vec (), chol_mat.rows (),
237 utmp.fortran_vec (), rw));
238 }
239 else
240 (*current_liboctave_error_handler) ("cholupdate: dimension mismatch"); 226 (*current_liboctave_error_handler) ("cholupdate: dimension mismatch");
227
228 FloatComplexColumnVector utmp = u;
229
230 OCTAVE_LOCAL_BUFFER (float, rw, n);
231
232 F77_XFCN (cch1up, CCH1UP, (n, chol_mat.fortran_vec (), chol_mat.rows (),
233 utmp.fortran_vec (), rw));
241 } 234 }
242 235
243 octave_idx_type 236 octave_idx_type
244 FloatComplexCHOL::downdate (const FloatComplexColumnVector& u) 237 FloatComplexCHOL::downdate (const FloatComplexColumnVector& u)
245 { 238 {
246 octave_idx_type info = -1; 239 octave_idx_type info = -1;
247 240
248 octave_idx_type n = chol_mat.rows (); 241 octave_idx_type n = chol_mat.rows ();
249 242
250 if (u.numel () == n) 243 if (u.numel () != n)
251 {
252 FloatComplexColumnVector utmp = u;
253
254 OCTAVE_LOCAL_BUFFER (float, rw, n);
255
256 F77_XFCN (cch1dn, CCH1DN, (n, chol_mat.fortran_vec (), chol_mat.rows (),
257 utmp.fortran_vec (), rw, info));
258 }
259 else
260 (*current_liboctave_error_handler) ("cholupdate: dimension mismatch"); 244 (*current_liboctave_error_handler) ("cholupdate: dimension mismatch");
245
246 FloatComplexColumnVector utmp = u;
247
248 OCTAVE_LOCAL_BUFFER (float, rw, n);
249
250 F77_XFCN (cch1dn, CCH1DN, (n, chol_mat.fortran_vec (), chol_mat.rows (),
251 utmp.fortran_vec (), rw, info));
261 252
262 return info; 253 return info;
263 } 254 }
264 255
265 octave_idx_type 256 octave_idx_type
270 261
271 octave_idx_type n = chol_mat.rows (); 262 octave_idx_type n = chol_mat.rows ();
272 263
273 if (u.numel () != n + 1) 264 if (u.numel () != n + 1)
274 (*current_liboctave_error_handler) ("cholinsert: dimension mismatch"); 265 (*current_liboctave_error_handler) ("cholinsert: dimension mismatch");
275 else if (j < 0 || j > n) 266 if (j < 0 || j > n)
276 (*current_liboctave_error_handler) ("cholinsert: index out of range"); 267 (*current_liboctave_error_handler) ("cholinsert: index out of range");
277 else 268
278 { 269 FloatComplexColumnVector utmp = u;
279 FloatComplexColumnVector utmp = u; 270
280 271 OCTAVE_LOCAL_BUFFER (float, rw, n);
281 OCTAVE_LOCAL_BUFFER (float, rw, n); 272
282 273 chol_mat.resize (n+1, n+1);
283 chol_mat.resize (n+1, n+1); 274
284 275 F77_XFCN (cchinx, CCHINX, (n, chol_mat.fortran_vec (), chol_mat.rows (),
285 F77_XFCN (cchinx, CCHINX, (n, chol_mat.fortran_vec (), chol_mat.rows (), 276 j + 1, utmp.fortran_vec (), rw, info));
286 j + 1, utmp.fortran_vec (), rw, info));
287 }
288 277
289 return info; 278 return info;
290 } 279 }
291 280
292 void 281 void
294 { 283 {
295 octave_idx_type n = chol_mat.rows (); 284 octave_idx_type n = chol_mat.rows ();
296 285
297 if (j < 0 || j > n-1) 286 if (j < 0 || j > n-1)
298 (*current_liboctave_error_handler) ("choldelete: index out of range"); 287 (*current_liboctave_error_handler) ("choldelete: index out of range");
299 else 288
300 { 289 OCTAVE_LOCAL_BUFFER (float, rw, n);
301 OCTAVE_LOCAL_BUFFER (float, rw, n); 290
302 291 F77_XFCN (cchdex, CCHDEX, (n, chol_mat.fortran_vec (), chol_mat.rows (),
303 F77_XFCN (cchdex, CCHDEX, (n, chol_mat.fortran_vec (), chol_mat.rows (), 292 j + 1, rw));
304 j + 1, rw)); 293
305 294 chol_mat.resize (n-1, n-1);
306 chol_mat.resize (n-1, n-1);
307 }
308 } 295 }
309 296
310 void 297 void
311 FloatComplexCHOL::shift_sym (octave_idx_type i, octave_idx_type j) 298 FloatComplexCHOL::shift_sym (octave_idx_type i, octave_idx_type j)
312 { 299 {
313 octave_idx_type n = chol_mat.rows (); 300 octave_idx_type n = chol_mat.rows ();
314 301
315 if (i < 0 || i > n-1 || j < 0 || j > n-1) 302 if (i < 0 || i > n-1 || j < 0 || j > n-1)
316 (*current_liboctave_error_handler) ("cholshift: index out of range"); 303 (*current_liboctave_error_handler) ("cholshift: index out of range");
317 else 304
318 { 305 OCTAVE_LOCAL_BUFFER (FloatComplex, w, n);
319 OCTAVE_LOCAL_BUFFER (FloatComplex, w, n); 306 OCTAVE_LOCAL_BUFFER (float, rw, n);
320 OCTAVE_LOCAL_BUFFER (float, rw, n); 307
321 308 F77_XFCN (cchshx, CCHSHX, (n, chol_mat.fortran_vec (), chol_mat.rows (),
322 F77_XFCN (cchshx, CCHSHX, (n, chol_mat.fortran_vec (), chol_mat.rows (), 309 i + 1, j + 1, w, rw));
323 i + 1, j + 1, w, rw));
324 }
325 } 310 }
326 311
327 #else 312 #else
328 313
329 void 314 void
331 { 316 {
332 warn_qrupdate_once (); 317 warn_qrupdate_once ();
333 318
334 octave_idx_type n = chol_mat.rows (); 319 octave_idx_type n = chol_mat.rows ();
335 320
336 if (u.length () == n) 321 if (u.length () != n)
337 {
338 init (chol_mat.hermitian () * chol_mat
339 + FloatComplexMatrix (u) * FloatComplexMatrix (u).hermitian (),
340 true, false);
341 }
342 else
343 (*current_liboctave_error_handler) ("cholupdate: dimension mismatch"); 322 (*current_liboctave_error_handler) ("cholupdate: dimension mismatch");
323
324 init (chol_mat.hermitian () * chol_mat
325 + FloatComplexMatrix (u) * FloatComplexMatrix (u).hermitian (),
326 true, false);
344 } 327 }
345 328
346 static bool 329 static bool
347 singular (const FloatComplexMatrix& a) 330 singular (const FloatComplexMatrix& a)
348 { 331 {
358 341
359 octave_idx_type info = -1; 342 octave_idx_type info = -1;
360 343
361 octave_idx_type n = chol_mat.rows (); 344 octave_idx_type n = chol_mat.rows ();
362 345
363 if (u.length () == n) 346 if (u.length () != n)
364 { 347 (*current_liboctave_error_handler) ("cholupdate: dimension mismatch");
365 if (singular (chol_mat)) 348
366 info = 2; 349 if (singular (chol_mat))
367 else 350 info = 2;
368 {
369 info = init (chol_mat.hermitian () * chol_mat
370 - FloatComplexMatrix (u)
371 * FloatComplexMatrix (u).hermitian (),
372 true, false);
373 if (info) info = 1;
374 }
375 }
376 else 351 else
377 (*current_liboctave_error_handler) ("cholupdate: dimension mismatch"); 352 {
353 info = init (chol_mat.hermitian () * chol_mat
354 - FloatComplexMatrix (u)
355 * FloatComplexMatrix (u).hermitian (),
356 true, false);
357 if (info) info = 1;
358 }
378 359
379 return info; 360 return info;
380 } 361 }
381 362
382 octave_idx_type 363 octave_idx_type
389 370
390 octave_idx_type n = chol_mat.rows (); 371 octave_idx_type n = chol_mat.rows ();
391 372
392 if (u.length () != n + 1) 373 if (u.length () != n + 1)
393 (*current_liboctave_error_handler) ("cholinsert: dimension mismatch"); 374 (*current_liboctave_error_handler) ("cholinsert: dimension mismatch");
394 else if (j < 0 || j > n) 375 if (j < 0 || j > n)
395 (*current_liboctave_error_handler) ("cholinsert: index out of range"); 376 (*current_liboctave_error_handler) ("cholinsert: index out of range");
377
378 if (singular (chol_mat))
379 info = 2;
380 else if (u(j).imag () != 0.0)
381 info = 3;
396 else 382 else
397 { 383 {
398 if (singular (chol_mat)) 384 FloatComplexMatrix a = chol_mat.hermitian () * chol_mat;
399 info = 2; 385 FloatComplexMatrix a1 (n+1, n+1);
400 else if (u(j).imag () != 0.0) 386 for (octave_idx_type k = 0; k < n+1; k++)
401 info = 3; 387 for (octave_idx_type l = 0; l < n+1; l++)
402 else 388 {
403 { 389 if (l == j)
404 FloatComplexMatrix a = chol_mat.hermitian () * chol_mat; 390 a1(k, l) = u(k);
405 FloatComplexMatrix a1 (n+1, n+1); 391 else if (k == j)
406 for (octave_idx_type k = 0; k < n+1; k++) 392 a1(k, l) = std::conj (u(l));
407 for (octave_idx_type l = 0; l < n+1; l++) 393 else
408 { 394 a1(k, l) = a(k < j ? k : k-1, l < j ? l : l-1);
409 if (l == j) 395 }
410 a1(k, l) = u(k); 396 info = init (a1, true, false);
411 else if (k == j) 397 if (info) info = 1;
412 a1(k, l) = std::conj (u(l));
413 else
414 a1(k, l) = a(k < j ? k : k-1, l < j ? l : l-1);
415 }
416 info = init (a1, true, false);
417 if (info) info = 1;
418 }
419 } 398 }
420 399
421 return info; 400 return info;
422 } 401 }
423 402
428 407
429 octave_idx_type n = chol_mat.rows (); 408 octave_idx_type n = chol_mat.rows ();
430 409
431 if (j < 0 || j > n-1) 410 if (j < 0 || j > n-1)
432 (*current_liboctave_error_handler) ("choldelete: index out of range"); 411 (*current_liboctave_error_handler) ("choldelete: index out of range");
433 else 412
434 { 413 FloatComplexMatrix a = chol_mat.hermitian () * chol_mat;
435 FloatComplexMatrix a = chol_mat.hermitian () * chol_mat; 414 a.delete_elements (1, idx_vector (j));
436 a.delete_elements (1, idx_vector (j)); 415 a.delete_elements (0, idx_vector (j));
437 a.delete_elements (0, idx_vector (j)); 416 init (a, true, false);
438 init (a, true, false);
439 }
440 } 417 }
441 418
442 void 419 void
443 FloatComplexCHOL::shift_sym (octave_idx_type i, octave_idx_type j) 420 FloatComplexCHOL::shift_sym (octave_idx_type i, octave_idx_type j)
444 { 421 {
446 423
447 octave_idx_type n = chol_mat.rows (); 424 octave_idx_type n = chol_mat.rows ();
448 425
449 if (i < 0 || i > n-1 || j < 0 || j > n-1) 426 if (i < 0 || i > n-1 || j < 0 || j > n-1)
450 (*current_liboctave_error_handler) ("cholshift: index out of range"); 427 (*current_liboctave_error_handler) ("cholshift: index out of range");
451 else 428
452 { 429 FloatComplexMatrix a = chol_mat.hermitian () * chol_mat;
453 FloatComplexMatrix a = chol_mat.hermitian () * chol_mat; 430 Array<octave_idx_type> p (dim_vector (n, 1));
454 Array<octave_idx_type> p (dim_vector (n, 1)); 431 for (octave_idx_type k = 0; k < n; k++) p(k) = k;
455 for (octave_idx_type k = 0; k < n; k++) p(k) = k; 432 if (i < j)
456 if (i < j) 433 {
457 { 434 for (octave_idx_type k = i; k < j; k++) p(k) = k+1;
458 for (octave_idx_type k = i; k < j; k++) p(k) = k+1; 435 p(j) = i;
459 p(j) = i; 436 }
460 } 437 else if (j < i)
461 else if (j < i) 438 {
462 { 439 p(j) = i;
463 p(j) = i; 440 for (octave_idx_type k = j+1; k < i+1; k++) p(k) = k-1;
464 for (octave_idx_type k = j+1; k < i+1; k++) p(k) = k-1; 441 }
465 } 442
466 443 init (a.index (idx_vector (p), idx_vector (p)), true, false);
467 init (a.index (idx_vector (p), idx_vector (p)), true, false);
468 }
469 } 444 }
470 445
471 #endif 446 #endif
472 447
473 FloatComplexMatrix 448 FloatComplexMatrix