Mercurial > hg > octave-nkf

comparison liboctave/sparse-base-chol.cc @ 7637:2be056f03720

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Fix fall back from sparse cholesky factorization to LU when matrix detected as not being positive definite
author David Bateman <dbateman@free.fr>
date Tue, 25 Mar 2008 16:17:16 -0400
parents daff886a8e2a
children 2176f2b4599e 72830070a17b
comparison
equal deleted inserted replaced
7636:99c410f7f0b0 7637:2be056f03720
163 163
164 cholmod_factor *Lfactor; 164 cholmod_factor *Lfactor;
165 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 165 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
166 Lfactor = CHOLMOD_NAME(analyze) (ac, cm); 166 Lfactor = CHOLMOD_NAME(analyze) (ac, cm);
167 CHOLMOD_NAME(factorize) (ac, Lfactor, cm); 167 CHOLMOD_NAME(factorize) (ac, Lfactor, cm);
168 cond = CHOLMOD_NAME(rcond) (Lfactor, cm);
169 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 168 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
170 169
171 minor_p = Lfactor->minor;
172 is_pd = cm->status == CHOLMOD_OK; 170 is_pd = cm->status == CHOLMOD_OK;
173 info = (is_pd ? 0 : cm->status); 171 info = (is_pd ? 0 : cm->status);
174 172
175 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 173 if (is_pd)
176 Lsparse = CHOLMOD_NAME(factor_to_sparse) (Lfactor, cm); 174 {
177 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
178
179 if (minor_p > 0 && minor_p < a_nr)
180 {
181 size_t n1 = a_nr + 1;
182 Lsparse->p = CHOLMOD_NAME(realloc) (minor_p+1, sizeof(octave_idx_type),
183 Lsparse->p, &n1, cm);
184 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 175 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
185 CHOLMOD_NAME(reallocate_sparse) 176 cond = CHOLMOD_NAME(rcond) (Lfactor, cm);
186 (static_cast<octave_idx_type *>(Lsparse->p)[minor_p], Lsparse, cm);
187 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 177 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
188 Lsparse->ncol = minor_p; 178
189 } 179 minor_p = Lfactor->minor;
190 180
191 drop_zeros (Lsparse); 181 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
192 182 Lsparse = CHOLMOD_NAME(factor_to_sparse) (Lfactor, cm);
193 if (! natural) 183 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
194 { 184
195 perms.resize (a_nr); 185 if (minor_p > 0 && minor_p < a_nr)
196 for (octave_idx_type i = 0; i < a_nr; i++) 186 {
197 perms(i) = static_cast<octave_idx_type *>(Lfactor->Perm)[i]; 187 size_t n1 = a_nr + 1;
198 } 188 Lsparse->p = CHOLMOD_NAME(realloc) (minor_p+1,
199 189 sizeof(octave_idx_type),
200 static char tmp[] = " "; 190 Lsparse->p, &n1, cm);
201 191 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
202 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 192 CHOLMOD_NAME(reallocate_sparse)
203 CHOLMOD_NAME(free_factor) (&Lfactor, cm); 193 (static_cast<octave_idx_type *>(Lsparse->p)[minor_p], Lsparse, cm);
204 CHOLMOD_NAME(finish) (cm); 194 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
205 CHOLMOD_NAME(print_common) (tmp, cm); 195 Lsparse->ncol = minor_p;
206 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 196 }
197
198 drop_zeros (Lsparse);
199
200 if (! natural)
201 {
202 perms.resize (a_nr);
203 for (octave_idx_type i = 0; i < a_nr; i++)
204 perms(i) = static_cast<octave_idx_type *>(Lfactor->Perm)[i];
205 }
206
207 static char tmp[] = " ";
208
209 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
210 CHOLMOD_NAME(free_factor) (&Lfactor, cm);
211 CHOLMOD_NAME(finish) (cm);
212 CHOLMOD_NAME(print_common) (tmp, cm);
213 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
214 }
207 #else 215 #else
208 (*current_liboctave_error_handler) 216 (*current_liboctave_error_handler)
209 ("Missing CHOLMOD. Sparse cholesky factorization disabled"); 217 ("Missing CHOLMOD. Sparse cholesky factorization disabled");
210 #endif 218 #endif
211 return info; 219 return info;