Hi again,

>>row (U, 0) = f;

> Aargh, I think the brand of coffee we are drinking, doen't work as expected :-)

Maybe you should try to put some cognac or whisky in your coffee ;o)

>>But garbage in U can't be a reason for getrs() to die.

>>How does it die? Error message? Program hangs?

>>Program aborts? Or you simply get garbage?

> It just aborts! For some weird reason, we cannot get VC71 to enter debug mode

> anymore. It says that the generated symbol file is corrupt. And we have no clue

> why?

When you had `row (U, 1) = f;' part of the memory which maybe

contained program code was overwritten.

>>Why do you invert stiffness matrix?

>>Numerical analysis books recommend back substitution:

>>It is more numerically stable and probably faster.

> It gave a tremendous speed-up. From 18 seconds downto 5 seconds.

Your stiffness matrix is probably symmetric positive definite.

So you can try potrf()/potrs().

> Thanks.

You're welcome ;o)

> Now, we "just" need to drop atlas and go for a sparse solution instead.

>

> What is the big difference between compressed_matrix and sparse_matrix anyhow?

sparse_matrix uses some kind of associative array to store

elements (e.g. ublas::map_array or std::map).

compressed_matrix uses three arrays -- for row major storage

`pointers' (sorted) to first non-zero element in each row,

column indices (sorted) of non-zero elements in each row,

non-zero elements.

sparse_matrix with std::map allows O(log nnz) random insertations

(where nnz is number of non-zero elements). Insertations in other

matrix types require moving all elements which are after the insert

position. So, for e.g. global stiffness matrix assembly sparse_matrix

with std::map is `ideal' data structure -- assembly is O(nnz * log nnz)

operation. (There is also coordinate_matrix. It also uses three arrays:

row and column indices of non-zero elements and non-zero elements.

When elements are inserted, they are in fact appended at the end (O(1)

operation), so matrix become unsorted. But before other uses, matrix

must be sorted -- again O(nnz * log nnz) operation.)

On the other hand, matrix traversals are usually faster for

compressed matrix -- storage is contiguous and there's no

need to `chase' pointer. And, probably even more important,

Fortran and C libraries for solving linear systems with sparse

matrices (e.g. SuperLU) usually expect compressed (or, sometimes,

coordinate) storage formats.

So, you can use sparse_matrix<std::map> or coordinate_matrix

for stiffness matrix assembly and then copy it to compressed_matrix

to solve the system.

You can find more information about sparse matrices in excellent

book by Yousef Saad `Iterative methods for sparse linear systems'

-- free .ps and .pdf can be found at

http://www-users.cs.umn.edu/~saad/
(Unfortunately, as Saad uses Fortran, associative array approach

is not mentioned -- not that associative arrays are impossible in

Fortran, but it's not simple to implement them; and there are

other, equally efficient and more `natural' solutions.)

Regards,

fres