lists@vrbk... lists@vrbk...
Tue Apr 21 16:42:01 CDT 2009

```hi guys,

in my program for solving the integral equations in the liquid state theory
(based on scipy/numpy; check pyoz.vrbka.net if you are interested in
details) i need to implement a solver of a linear system AX=B - generally,
this shouldn't be a big problem. there are some issues, however...

1) X and B are sets of functions - X_ij and B_ij, X and B are NxN matrices
(most frequently N=2-3)
2) the functions are discretized (generally 2^10-2^12 discretization
points)

therefore, i have (in the easiest case) 2^10 of 2x2 problems

3) A is slightly complicated linear operator, involving fourier transforms

in the literature, i found that this system of equations could be solved
using the nonsymmetric (iterative) conjugate gradients procedure. in the
paper (i can provide the reference on request, don't have it here at the
moment), actually, the author presented the equations for the CG algorithm,
involving the adjoint of the operator A - the resulting equations were
relatively simple. the author speaks also of matrices of the size of 10^3 x
10^3 without giving further details - i presume functions at all the
discretization points would be put into one big block diagonal matrix...

having a look at the sparse.linalg.solve.cg, there's no way of specifying
the adjoint (which should be, if i understand the method correctly,
necessary for getting the result, since it's impossible to invert something
that big efficiently).

i'm searching for help here, since my knowledge of the mathematics is quite
limited in this area... would be really very glad for any hints or comments
- how would you proceed in this case?