On 25 May 2010, at 21:01, Jaroslav Hajek wrote: > On Tue, May 25, 2010 at 6:00 PM, Carlo de Falco <[email protected] > > wrote: >> >> On 25 May 2010, at 16:17, Jaroslav Hajek wrote: >> >>> On Tue, May 25, 2010 at 3:13 PM, Carlo de Falco <[email protected] >>> > >>> wrote: >>>> >>>> On 25 May 2010, at 14:59, Jaroslav Hajek wrote: >>>> >>>>> unless anyone objects, I'll become the maintainer of this package. >>>> >>>> oh, well.. as now linear algebra has a maintainer I finally have >>>> someone >>>> to >>>> direct my complaints to ;) >>>> >>>>> comments, further suggestions? >>>> >>>> well yes, some time ago an initial implementation of pgmres was >>>> added to >>>> the >>>> package but the author[*] >>>> seems to not have taken too much care for making it compatibile >>>> with >>>> matlab's gmres function, maybe you >>>> could have a look into this before the release? >>>> >>>> c. >>>> >>>> [*] OK, you got me, that lazy author is me... pgmres needs a lot of >>>> improvement but it has been on my todo list >>>> for a long time now, so although I think it would make sense to >>>> complete >>>> it >>>> before the release I will not at all >>>> ask to delay the release because of this. >>>> >>> >>> >>> Is there particular reason why pgmres is a C++ file? >> >> the actual reason why pgmres was written in the first place was to >> demonstrate the overloading of the '*' >> operator in c++ and the use of a function instead of a matrix in the >> matrix-free implementation of an iterative solver. >> In principle I have no objections to switching to an m file if >> performance >> is not affected, though. >> >>> The only inner >>> loop I can see is the modified Gram-Schmidt orthogonalization: >>> >>> H.resize (reit+1, reit, 0.0); >>> for (octave_idx_type j = 0; j < reit; j++) >>> { >>> H(j, reit-1) = (V.column (j).transpose ()) * tmp; >>> tmp = tmp - (H (j, reit-1) * V.column (j)); >>> } >>> >>> which would be useful to have as a separate function, such as: >>> >>> >>> function [y, r] = mgorth (v, x); >>> # orthogonalizes a column vector x relative to orthogonal basis v >>> using the modified gram-schmidt orthogonalization. >>> # mathematically, this is equivalent to >>> # r = v'*x; y = x - v*r; >>> # but computed in a more stable manner that avoids accumulating >>> roundoff errors when done sequentially. >>> >>> >>> pgmres could then be written as a m-file. What do you think? >> >> why is only the inner loop important, if many iterations are >> performed >> doesn't the outer loop also impact speed? >> > > If a loop contains "heavy" operations, such as large matrix > multiplication, matrix division, factorization, convolution, fft, > etc., the interpreter overhead will typically become negligible > compared to the time taken by these operations, so it doesn't matter > whether it's written in m-code or C++. > Also, when a callback function is called in the loop, it will most > likely be an m-function, so the loop will involve interpreter overhead > anyway.
how would you do the matrix-free matrix matrix vector multiplication in an m-file? would if (isa (A, "function_handle") || ischar (A)) y = feval (A, x); else y = A * x; endif be as efficient as the overloaded "*"? As I said, I initially wrote pgmres as a teaching example that needed to be in C++ but if its efficiency is not much reduced, I agree making it an m-file as that would make it easier to maintain. I'll try to do this c. ------------------------------------------------------------------------------ _______________________________________________ Octave-dev mailing list [email protected] https://lists.sourceforge.net/lists/listinfo/octave-dev
