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-unsubscr...@commons.apache.org For additional commands, e-mail: dev-h...@commons.apache.org
