Mercurial > octave-libtiff

comparison liboctave/numeric/sparse-chol.cc @ 21177:a10f60e13243

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
style fixes in liboctave/numeric directory * base-qr.h, randgamma.h, randmtzig.h, randpoisson.h, sparse-chol.cc, sparse-chol.h, sparse-dmsolve.cc, sparse-lu.cc, sparse-lu.h, sparse-qr.cc, sparse-qr.h, oct-sparse.h: Style fixes.
author John W. Eaton <jwe@octave.org>
date Tue, 02 Feb 2016 16:53:01 -0500
parents 307096fb67e1
children 7f35125714b4
comparison
equal deleted inserted replaced
21176:791dcb32b657 21177:a10f60e13243
1 /* 1 /*
2 2
3 Copyright (C) 2016 John W. Eaton
3 Copyright (C) 2005-2015 David Bateman 4 Copyright (C) 2005-2015 David Bateman
4 Copyright (C) 1998-2005 Andy Adler 5 Copyright (C) 1998-2005 Andy Adler
5 6
6 This file is part of Octave. 7 This file is part of Octave.
7 8
135 136
136 template <typename chol_type> 137 template <typename chol_type>
137 void 138 void
138 sparse_chol<chol_type>::sparse_chol_rep::drop_zeros (const cholmod_sparse *S) 139 sparse_chol<chol_type>::sparse_chol_rep::drop_zeros (const cholmod_sparse *S)
139 { 140 {
140 chol_elt sik;
141 octave_idx_type *Sp, *Si;
142 chol_elt *Sx;
143 octave_idx_type pdest, k, ncol, p, pend;
144
145 if (! S) 141 if (! S)
146 return; 142 return;
147 143
148 Sp = static_cast<octave_idx_type *>(S->p); 144 octave_idx_type *Sp = static_cast<octave_idx_type *>(S->p);
149 Si = static_cast<octave_idx_type *>(S->i); 145 octave_idx_type *Si = static_cast<octave_idx_type *>(S->i);
150 Sx = static_cast<chol_elt *>(S->x); 146 chol_elt *Sx = static_cast<chol_elt *>(S->x);
151 pdest = 0; 147
152 ncol = S->ncol; 148 octave_idx_type pdest = 0;
153 149 octave_idx_type ncol = S->ncol;
154 for (k = 0; k < ncol; k++) 150
155 { 151 for (octave_idx_type k = 0; k < ncol; k++)
156 p = Sp[k]; 152 {
157 pend = Sp[k+1]; 153 octave_idx_type p = Sp[k];
154 octave_idx_type pend = Sp[k+1];
158 Sp[k] = pdest; 155 Sp[k] = pdest;
156
159 for (; p < pend; p++) 157 for (; p < pend; p++)
160 { 158 {
161 sik = Sx[p]; 159 chol_elt sik = Sx[p];
160
162 if (CHOLMOD_IS_NONZERO (sik)) 161 if (CHOLMOD_IS_NONZERO (sik))
163 { 162 {
164 if (p != pdest) 163 if (p != pdest)
165 { 164 {
166 Si[pdest] = Si[p]; 165 Si[pdest] = Si[p];
167 Sx[pdest] = sik; 166 Sx[pdest] = sik;
168 } 167 }
168
169 pdest++; 169 pdest++;
170 } 170 }
171 } 171 }
172 } 172 }
173
173 Sp[ncol] = pdest; 174 Sp[ncol] = pdest;
174 } 175 }
175 176
176 // Must provide a specialization for this function. 177 // Must provide a specialization for this function.
177 template <typename T> 178 template <typename T>
200 bool natural, bool force) 201 bool natural, bool force)
201 { 202 {
202 volatile octave_idx_type info = 0; 203 volatile octave_idx_type info = 0;
203 204
204 #ifdef HAVE_CHOLMOD 205 #ifdef HAVE_CHOLMOD
206
205 octave_idx_type a_nr = a.rows (); 207 octave_idx_type a_nr = a.rows ();
206 octave_idx_type a_nc = a.cols (); 208 octave_idx_type a_nc = a.cols ();
207 209
208 if (a_nr != a_nc) 210 if (a_nr != a_nc)
209 (*current_liboctave_error_handler) ("sparse_chol requires square matrix"); 211 (*current_liboctave_error_handler) ("sparse_chol requires square matrix");
210 212
211 cholmod_common *cm = &Common; 213 cholmod_common *cm = &Common;
212 214
213 // Setup initial parameters 215 // Setup initial parameters
216
214 CHOLMOD_NAME(start) (cm); 217 CHOLMOD_NAME(start) (cm);
215 cm->prefer_zomplex = false; 218 cm->prefer_zomplex = false;
216 219
217 double spu = octave_sparse_params::get_key ("spumoni"); 220 double spu = octave_sparse_params::get_key ("spumoni");
221
218 if (spu == 0.) 222 if (spu == 0.)
219 { 223 {
220 cm->print = -1; 224 cm->print = -1;
221 SUITESPARSE_ASSIGN_FPTR (printf_func, cm->print_function, 0); 225 SUITESPARSE_ASSIGN_FPTR (printf_func, cm->print_function, 0);
222 } 226 }
225 cm->print = static_cast<int> (spu) + 2; 229 cm->print = static_cast<int> (spu) + 2;
226 SUITESPARSE_ASSIGN_FPTR (printf_func, cm->print_function, &SparseCholPrint); 230 SUITESPARSE_ASSIGN_FPTR (printf_func, cm->print_function, &SparseCholPrint);
227 } 231 }
228 232
229 cm->error_handler = &SparseCholError; 233 cm->error_handler = &SparseCholError;
234
230 SUITESPARSE_ASSIGN_FPTR2 (divcomplex_func, cm->complex_divide, divcomplex); 235 SUITESPARSE_ASSIGN_FPTR2 (divcomplex_func, cm->complex_divide, divcomplex);
231 SUITESPARSE_ASSIGN_FPTR2 (hypot_func, cm->hypotenuse, hypot); 236 SUITESPARSE_ASSIGN_FPTR2 (hypot_func, cm->hypotenuse, hypot);
232 237
233 cm->final_asis = false; 238 cm->final_asis = false;
234 cm->final_super = false; 239 cm->final_super = false;
238 cm->final_resymbol = false; 243 cm->final_resymbol = false;
239 244
240 cholmod_sparse A; 245 cholmod_sparse A;
241 cholmod_sparse *ac = &A; 246 cholmod_sparse *ac = &A;
242 double dummy; 247 double dummy;
248
243 ac->nrow = a_nr; 249 ac->nrow = a_nr;
244 ac->ncol = a_nc; 250 ac->ncol = a_nc;
245 251
246 ac->p = a.cidx (); 252 ac->p = a.cidx ();
247 ac->i = a.ridx (); 253 ac->i = a.ridx ();
300 Lsparse->p, &n1, cm); 306 Lsparse->p, &n1, cm);
301 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 307 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
302 CHOLMOD_NAME(reallocate_sparse) 308 CHOLMOD_NAME(reallocate_sparse)
303 (static_cast<octave_idx_type *>(Lsparse->p)[minor_p], Lsparse, cm); 309 (static_cast<octave_idx_type *>(Lsparse->p)[minor_p], Lsparse, cm);
304 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 310 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
311
305 Lsparse->ncol = minor_p; 312 Lsparse->ncol = minor_p;
306 } 313 }
307 314
308 drop_zeros (Lsparse); 315 drop_zeros (Lsparse);
309 316
334 template <typename chol_type> 341 template <typename chol_type>
335 SparseMatrix 342 SparseMatrix
336 sparse_chol<chol_type>::sparse_chol_rep::Q (void) const 343 sparse_chol<chol_type>::sparse_chol_rep::Q (void) const
337 { 344 {
338 #ifdef HAVE_CHOLMOD 345 #ifdef HAVE_CHOLMOD
346
339 octave_idx_type n = Lsparse->nrow; 347 octave_idx_type n = Lsparse->nrow;
340 SparseMatrix p (n, n, n); 348 SparseMatrix p (n, n, n);
341 349
342 for (octave_idx_type i = 0; i < n; i++) 350 for (octave_idx_type i = 0; i < n; i++)
343 { 351 {
344 p.xcidx (i) = i; 352 p.xcidx (i) = i;
345 p.xridx (i) = static_cast<octave_idx_type>(perms (i)); 353 p.xridx (i) = static_cast<octave_idx_type>(perms (i));
346 p.xdata (i) = 1; 354 p.xdata (i) = 1;
347 } 355 }
356
348 p.xcidx (n) = n; 357 p.xcidx (n) = n;
349 358
350 return p; 359 return p;
360
351 #else 361 #else
362
352 return SparseMatrix (); 363 return SparseMatrix ();
364
353 #endif 365 #endif
354 } 366 }
355 367
356 template <typename chol_type> 368 template <typename chol_type>
357 sparse_chol<chol_type>::sparse_chol (void) 369 sparse_chol<chol_type>::sparse_chol (void)
418 template <typename chol_type> 430 template <typename chol_type>
419 chol_type 431 chol_type
420 sparse_chol<chol_type>::L (void) const 432 sparse_chol<chol_type>::L (void) const
421 { 433 {
422 #ifdef HAVE_CHOLMOD 434 #ifdef HAVE_CHOLMOD
435
423 cholmod_sparse *m = rep->L (); 436 cholmod_sparse *m = rep->L ();
437
424 octave_idx_type nc = m->ncol; 438 octave_idx_type nc = m->ncol;
425 octave_idx_type nnz = m->nzmax; 439 octave_idx_type nnz = m->nzmax;
440
426 chol_type ret (m->nrow, nc, nnz); 441 chol_type ret (m->nrow, nc, nnz);
442
427 for (octave_idx_type j = 0; j < nc+1; j++) 443 for (octave_idx_type j = 0; j < nc+1; j++)
428 ret.xcidx (j) = static_cast<octave_idx_type *>(m->p)[j]; 444 ret.xcidx (j) = static_cast<octave_idx_type *>(m->p)[j];
445
429 for (octave_idx_type i = 0; i < nnz; i++) 446 for (octave_idx_type i = 0; i < nnz; i++)
430 { 447 {
431 ret.xridx (i) = static_cast<octave_idx_type *>(m->i)[i]; 448 ret.xridx (i) = static_cast<octave_idx_type *>(m->i)[i];
432 ret.xdata (i) = static_cast<chol_elt *>(m->x)[i]; 449 ret.xdata (i) = static_cast<chol_elt *>(m->x)[i];
433 } 450 }
451
434 return ret; 452 return ret;
453
435 #else 454 #else
455
436 return chol_type (); 456 return chol_type ();
457
437 #endif 458 #endif
438 } 459 }
439 460
440 template <typename chol_type> 461 template <typename chol_type>
441 octave_idx_type 462 octave_idx_type
475 template <typename chol_type> 496 template <typename chol_type>
476 chol_type 497 chol_type
477 sparse_chol<chol_type>::inverse (void) const 498 sparse_chol<chol_type>::inverse (void) const
478 { 499 {
479 chol_type retval; 500 chol_type retval;
480 #ifdef HAVE_CHOLMOD 501
502 #ifdef HAVE_CHOLMOD
503
481 cholmod_sparse *m = rep->L (); 504 cholmod_sparse *m = rep->L ();
482 octave_idx_type n = m->ncol; 505 octave_idx_type n = m->ncol;
483 ColumnVector perms = rep->perm (); 506 ColumnVector perms = rep->perm ();
484 double rcond2; 507 double rcond2;
485 octave_idx_type info; 508 octave_idx_type info;
487 chol_type linv = L ().hermitian ().inverse (mattype, info, rcond2, 1, 0); 510 chol_type linv = L ().hermitian ().inverse (mattype, info, rcond2, 1, 0);
488 511
489 if (perms.numel () == n) 512 if (perms.numel () == n)
490 { 513 {
491 SparseMatrix Qc = Q (); 514 SparseMatrix Qc = Q ();
515
492 retval = Qc * linv * linv.hermitian () * Qc.transpose (); 516 retval = Qc * linv * linv.hermitian () * Qc.transpose ();
493 } 517 }
494 else 518 else
495 retval = linv * linv.hermitian (); 519 retval = linv * linv.hermitian ();
496 #endif 520
521 #endif
522
497 return retval; 523 return retval;
498 } 524 }
499 525
500 template <typename chol_type> 526 template <typename chol_type>
501 chol_type 527 chol_type