Mailing-List: contact dev-help@commons.apache.org; run by ezmlm
Precedence: bulk
Reply-To: "Commons Developers List" <dev@commons.apache.org>
Received-SPF: neutral (nike.apache.org: 129.104.30.34 is neither permitted nor
 denied by domain of sebastien.brisard@m4x.org)
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit
From: "=?UTF-8?Q?S=C3=A9bastien=20Brisard?=" <sebastien.brisard@m4x.org>
Subject: [math] Iterative linear solvers
To: dev@commons.apache.org
Message-Id: <20110413062230.2B0E41405982C@svoboda.polytechnique.org>
Date: Wed, 13 Apr 2011 08:22:30 +0200 (CEST)

Dear All,
first of all, I should say that I'm new to the Commons-Math community, and I
therefore apologize for any mistake/inconsistency.
I'm currently developing a numerical method for the simulation of highly
heterogeneous materials. This requires the inversion of a (big) linear system,
whose matrix A is *not* known in closed form. The only thing I know about A is
how to compute the product A.x for any vector x. True enough, this would allow
me to compute the coefficients of the matrix, but
1. I do not need them,
2. It would waste memory space.

Many iterative solvers actually do not require the coefficients of matrix A.
For example, I've implemented CG and SYMMLQ, only assuming that A is defined
by a method operate(double[] arg, double[] img).

If you feel that this code might find its place in Commons-Math, I'm very
willing to contribute it. However, here are a few questions
1. The matrix as I defined it does not fit into the
org.apache.commons.math.linear.RealMatrix interface (since I do not want to
define add, getRow, etc...). For the time being, I defined a so called
LinearMap interface, which inherits from AnyMatrix, but only defines one
additional method
void operate(double[] arg, double[] img)

2. BTW, the vector in which the image A.x of the argument x must be stored is
passed as an argument, in order to limit memory allocation. I always do that,
but noticed that it is not common-practice in Commons-Math. Should I give up
this habit? Is it bad-practice? Is the garbage collector so efficient nowadays
that such precaution has become pointless?

2bis. If you think I should substitute
double[] operate(double[] arg)
to
void operate(double[] arg, double[] img)
then RealMatrix might extend LinearMap (or whatever its name).

3. Regarding preconditionners: the GC or SYMMLQ solvers can take as an
argument a preconditionner M. Again, we do not need to know its coefficients.
The only thing which is required is to compute M.x or inv(M).y. So I have
defined another interface, namely InvertibleLinearMap, which inherits from
LinearMap, with one additional method
void solve(double[] b, double[] x)
Any thoughts?

3. Regarding potential exceptions that may be returned. I'm puzzled by the
constructors of the exceptions which might be of use for me, namely
NonSymmetricMatrixException
NonPositiveDefiniteMatrixException

both exceptions refer to row/column indices, which I do not know! For example,
I would like to return a NonPositiveDefiniteMatrixException if, during the
iterations, I find that transpose(x).A.x < 0, without referring to a specific
row or column. Besides, such an exception might be thrown either when M or A
is non-symmetric (both passed to the GC/SYMMLQ solver at the same time), and I
would like the message of the exception to state wether the problem came from
the matrix A itself, or the preconditionner M. How should I do? Should I
define my own exceptions? This seems to me to be a waste...

I'm looking forward to having the opportunity to sharing my code. Meanwhile, I
do apologize if the previous remarks and questions sound silly.

Best,
Sebastien

---------------------------------------------------------------------
To unsubscribe, e-mail: dev-unsubscribe@commons.apache.org
For additional commands, e-mail: dev-help@commons.apache.org