Le 13/04/2011 08:22, Sébastien Brisard a écrit : > 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
Yes, it would be interesting. > 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) This seems fine. As you wrote, computing all coefficients would be wasting resources, so just having a different simplified interface with is nice. > > 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? We sometimes provide an additional boolean argument for such cases, so the user can specify if the array should be cloned to be kept independent from the caller or if it can be referenced. > > 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). That seems a good idea. > > 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? No idea for now, we should have a look at it. > > 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 You can set up new constructors if needed with specific messages. Beware that in the current repository trunk (which heads towards 3.0), there have been many changes in the exception, so you'll have to make sure your code does fit with this new scheme. > 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. No worries! Also as we have a specific process to deal with incoming contributions, we will also need a Software Grant (see <http://www.apache.org/licenses/#grants>). Thanks Luc > > Best, > Sebastien > > --------------------------------------------------------------------- > To unsubscribe, e-mail: dev-unsubscr...@commons.apache.org > For additional commands, e-mail: dev-h...@commons.apache.org > > --------------------------------------------------------------------- To unsubscribe, e-mail: dev-unsubscr...@commons.apache.org For additional commands, e-mail: dev-h...@commons.apache.org
