On Tue, May 25, 2010 at 9:29 PM, Carlo de Falco <[email protected]> wrote:
>
> 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 "*"?
>
Or we could do it like this
if (ischar (A))
Ax = str2func (A);
elseif (ismatrix (A))
Ax = @(x) A*x;
elseif (isa (A, "function_handle")
Ax = A;
else
error ("...");
endif
and then do y = Ax (x). No doubt it will be slower than operator
overloading in C++, the question is whether it matters, given that the
operation wrapped is typically quite expensive. I have plans for
optimizing anonymous functions that only involve a simple expression,
but no code yet, so now it more or less involves an almost full
function call overhead (in fact it is slightly faster, but not much).
> 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.
OK. Please note that your C++ version only works in real doubles,
whereas a properly written m-file function will probably work well in
both single and double (that is typically automatic, one only needs to
be careful about eps() and tolerances) and both real and complex (that
is sometimes harder, requires using ' and .' properly, sometimes it is
not possible).
--
RNDr. Jaroslav Hajek, PhD
computing expert & GNU Octave developer
Aeronautical Research and Test Institute (VZLU)
Prague, Czech Republic
url: www.highegg.matfyz.cz
------------------------------------------------------------------------------
_______________________________________________
Octave-dev mailing list
[email protected]
https://lists.sourceforge.net/lists/listinfo/octave-dev
