[julia-users] Defining a type that inherits methods

2016-10-29 Thread Jérémy Béjanin
I know that it is not possible to subtype concrete types, but is it 
possible in some other way to create a type that behave exactly like an 
existing concrete type (in my case a Matrix) such that it would keep all 
the methods associated with it, but would still dispatch when a more 
specific (ie specific to that new type) method exists?

I searched on the mailing list but could not find anything.

Thanks,
Jeremy


Re: [julia-users] Parsing complex numbers

2016-10-28 Thread Jérémy Béjanin
How easy would it be to make one assuming the standard julia notation 3+4im?

On Friday, October 28, 2016 at 3:58:08 PM UTC-4, Yichao Yu wrote:
>
> On Fri, Oct 28, 2016 at 3:53 PM, Jérémy Béjanin 
> > wrote: 
> > I've noticed that parsing a string representing a real number yields a 
> real 
> > number, but parsing a string representing a complex number yields an 
> > expression that must subsequently be evaluated. Is there a reason for 
> that 
> > behaviour? I'd like to avoid that behaviour considering I am reading 
> > user-inputted data. 
> > 
> > ```julia 
> > julia> typeof(parse("1.60254+3im")) 
> > Expr 
> > 
> > julia> typeof(parse("1.60254")) 
> > Float64 
> > ``` 
>
> Do no use `parse(::String)` to parse numbers. It is for parsing 
> generic julia code. Use `parse(Float64, str)` to parse a floating 
> point number. 
>
> I don't think we have a parsing function for complex number likely 
> because there isn't a universal standard. 
>


[julia-users] Parsing complex numbers

2016-10-28 Thread Jérémy Béjanin
I've noticed that parsing a string representing a real number yields a real 
number, but parsing a string representing a complex number yields an 
expression that must subsequently be evaluated. Is there a reason for that 
behaviour? I'd like to avoid that behaviour considering I am reading 
user-inputted data.

```julia
julia> typeof(parse("1.60254+3im"))
Expr

julia> typeof(parse("1.60254"))
Float64
```


Re: [julia-users] matrix multiplications

2016-10-18 Thread Jérémy Béjanin
I know, I was asking about that being the default behaviour of *

On Tuesday, October 18, 2016 at 2:10:14 PM UTC-4, Stefan Karpinski wrote:
>
> That's what A_mul_B! does.
>
> On Tue, Oct 18, 2016 at 1:45 PM, Jérémy Béjanin  > wrote:
>
>> I think this is something I might have read about in the past, but are 
>> there plans to make y = a*b use an already allocated y?
>>
>> On Tuesday, October 18, 2016 at 12:38:00 PM UTC-4, Stefan Karpinski wrote:
>>>
>>>   A_mul_B!(Y, A, B) -> Y
>>>
>>>   Calculates the matrix-matrix or matrix-vector product A⋅B and stores 
>>> the result in Y, overwriting the
>>>   existing value of Y. Note that Y must not be aliased with either A or 
>>> B.
>>>
>>>   julia> A=[1.0 2.0; 3.0 4.0]; B=[1.0 1.0; 1.0 1.0]; Y = similar(B); 
>>> A_mul_B!(Y, A, B);
>>>
>>>   julia> Y
>>>   2×2 Array{Float64,2}:
>>>3.0  3.0
>>>7.0  7.0
>>>
>>> On Tue, Oct 18, 2016 at 10:27 AM,  wrote:
>>>
>>>> hi guys
>>>> is there a way to reduce allocated memory in matrix multiplications?
>>>>
>>>> for example this code blew in my machine gives :
>>>>
>>>> function test4(n)
>>>>   a = rand(n,n)
>>>>   for i = 1:100
>>>>  a*a
>>>>end
>>>>  end
>>>>
>>>>
>>>> -- answer  --
>>>> test4(1)
>>>> # force compiling 
>>>>
>>>> @time test4(1000)
>>>> 16.589743 seconds (433 allocations: 770.587 MB, 0.68% gc time)
>>>>
>>>
>>>
>

Re: [julia-users] Recasting arrays

2016-10-18 Thread Jérémy Béjanin
Thanks, not sure how I missed that...

On Tuesday, October 18, 2016 at 1:55:58 PM UTC-4, Yichao Yu wrote:
>
>
>
> On Tue, Oct 18, 2016 at 1:38 PM, Jérémy Béjanin  > wrote:
>
>> I have seen the rem(a,b) function being used to recast numbers, but I was 
>> wondering if there was a way to recast arrays as in C.
>>
>
> Do note that this is undefined in standard C.
>  
>
>>
>> Say, for example, that I have an array of bytes, is it possible to recast 
>> that array as UInt16, by taking pairs of bytes, or as Int32 by taking each 
>> consecutive 4 bytes?
>>
>
> reinterpret. Note the endianess dependency.
>  
>
>>
>> Thanks,
>> Jeremy
>>
>
>

Re: [julia-users] matrix multiplications

2016-10-18 Thread Jérémy Béjanin
I think this is something I might have read about in the past, but are 
there plans to make y = a*b use an already allocated y?

On Tuesday, October 18, 2016 at 12:38:00 PM UTC-4, Stefan Karpinski wrote:
>
>   A_mul_B!(Y, A, B) -> Y
>
>   Calculates the matrix-matrix or matrix-vector product A⋅B and stores the 
> result in Y, overwriting the
>   existing value of Y. Note that Y must not be aliased with either A or B.
>
>   julia> A=[1.0 2.0; 3.0 4.0]; B=[1.0 1.0; 1.0 1.0]; Y = similar(B); 
> A_mul_B!(Y, A, B);
>
>   julia> Y
>   2×2 Array{Float64,2}:
>3.0  3.0
>7.0  7.0
>
> On Tue, Oct 18, 2016 at 10:27 AM, > 
> wrote:
>
>> hi guys
>> is there a way to reduce allocated memory in matrix multiplications?
>>
>> for example this code blew in my machine gives :
>>
>> function test4(n)
>>   a = rand(n,n)
>>   for i = 1:100
>>  a*a
>>end
>>  end
>>
>>
>> -- answer  --
>> test4(1)
>> # force compiling 
>>
>> @time test4(1000)
>> 16.589743 seconds (433 allocations: 770.587 MB, 0.68% gc time)
>>
>
>

[julia-users] Recasting arrays

2016-10-18 Thread Jérémy Béjanin
I have seen the rem(a,b) function being used to recast numbers, but I was 
wondering if there was a way to recast arrays as in C.

Say, for example, that I have an array of bytes, is it possible to recast 
that array as UInt16, by taking pairs of bytes, or as Int32 by taking each 
consecutive 4 bytes?

Thanks,
Jeremy


Re: [julia-users] Calling DLL function

2016-10-07 Thread Jérémy Béjanin
The format is fairly complicated, but I need to take 4 bytes from that 
buffer and extract two 12-bit numbers from that 32-bit string (along with 
some reversing operations, as you can see in the function). I cannot seem 
to find bit operations that address/concatenate individual bits without 
going through a string format and parsing. Would you have any guidance as 
to where in the manual I should look?

On Thursday, October 6, 2016 at 8:15:41 PM UTC-4, Yichao Yu wrote:
>
> On Thu, Oct 6, 2016 at 5:40 PM, Jérémy Béjanin  > wrote: 
> > Thanks! this works perfectly. 
> > 
> > The pointer_to_array() function seems to be deprecated, however, do you 
> know 
> > how to use the suggested replacement, unsafe_wrap()? 
>
> Yes. 
>
> > 
> > And if it's not too much to ask, I am wondering how I can do the 
> conversion 
> > from this UInt8 buffer more efficiently. This is what I am currently 
> doing, 
> > where blocks is the buffer: 
> > 
> > # We now have all the blocks in the blocks array, we need to extract the 
> > # 12-bit samples from that array. The data in the blocks is organized in 
> 4 
> > # bytes / 32-bit words, and each word contains two 12-bit samples. 
> > # Calculate the number of 32-bit words in the blocks array 
> > numwords = div(numblocks*DIG_BLOCK_SIZE,4) 
> > # Initialize array to store the samples 
> > data = Array{Int32}(2*numwords) 
> > for n=1:numwords 
> > # Convert next four bytes to 32-bit string 
> > s = 
> > 
> reverse(bits(blocks[(n-1)*4+1]))*reverse(bits(blocks[(n-1)*4+2]))*reverse(bits(blocks[(n-1)*4+3]))*reverse(bits(blocks[(n-1)*4+4]))
>  
>
> > # Parse the first 12 bits as the first sample and the 12 bits 4 bits 
> > later as the second sample 
> > data[(n-1)*2+1] = parse(Int,reverse(s[1:12]),2) 
> > data[(n-1)*2+2] = parse(Int,reverse(s[17:28]),2) 
> > end 
> > 
> > This is pretty slow, I assume due to the translation between numbers and 
> > strings. Is there a better way to do this? 
>
> I don't know what exactly the format is but you should not go though 
> string and you should be able to do what you want with simple 
> interger/bits operations. 
>
> > 
> > Thanks, 
> > Jeremy 
> > 
>


Re: [julia-users] Calling DLL function

2016-10-06 Thread Jérémy Béjanin
Thanks! this works perfectly.

The pointer_to_array() function seems to be deprecated, however, do you 
know how to use the suggested replacement, unsafe_wrap()?

And if it's not too much to ask, I am wondering how I can do the conversion 
from this UInt8 buffer more efficiently. This is what I am currently doing, 
where blocks is the buffer:

# We now have all the blocks in the blocks array, we need to extract the
# 12-bit samples from that array. The data in the blocks is organized in 4
# bytes / 32-bit words, and each word contains two 12-bit samples.
# Calculate the number of 32-bit words in the blocks array
numwords = div(numblocks*DIG_BLOCK_SIZE,4)
# Initialize array to store the samples
data = Array{Int32}(2*numwords)
for n=1:numwords
# Convert next four bytes to 32-bit string
s = 
reverse(bits(blocks[(n-1)*4+1]))*reverse(bits(blocks[(n-1)*4+2]))*reverse(bits(blocks[(n-1)*4+3]))*reverse(bits(blocks[(n-1)*4+4]))
# Parse the first 12 bits as the first sample and the 12 bits 4 bits 
later as the second sample
data[(n-1)*2+1] = parse(Int,reverse(s[1:12]),2)
data[(n-1)*2+2] = parse(Int,reverse(s[17:28]),2)
end

This is pretty slow, I assume due to the translation between numbers and 
strings. Is there a better way to do this?

Thanks,
Jeremy



[julia-users] Calling DLL function

2016-10-05 Thread Jérémy Béjanin
Hello,

I am trying to call a DLL function to allocate a page aligned array from 
julia. Here is a C++ example using that function. Note that the DLL 
function is wrapped into this uvAPI api, but is otherwise the same except 
for the name. The DLL function is called x_MemAlloc, while the API function 
is called X_MemAlloc. The function should return 0 on success, or 1 on 
failure.

#define DIG_BLOCK_SIZE 1024 * 1024

uvAPI 

 
*uv = new uvAPI 

;
unsigned char * sysMem = NULL;
error = uv->X_MemAlloc 

((void**)&sysMem, DIG_BLOCK_SIZE 

);
if (error)
{
 std::cout << "failed to allocate block buffer" << std::endl;
 return  -1;
}

Here is my julia code:

DIG_BLOCK_SIZE = 1024 * 1024
block = Array{Cuchar}(DIG_BLOCK_SIZE)
if ccall((:x_MemAlloc,AcqSynth),Cint,(Ptr{Ptr{Void}},Csize_t),pointer(block
),DIG_BLOCK_SIZE)==1
error("Failed to allocate block buffer!")
end

However, the function returns 1 whether I use Ptr{Void} and block, or 
Ptr{Ptr{Void}} and pointer(block). Is this the proper translation?

I am trying the acquire data with a PCIe digitizer and currently reading 
out the data into the buffer does not work. I am wondering if this is 
because the buffer is not allocated properly.

Note that other DLL functions I've tried work properly.

Thanks,
Jeremy


Re: [julia-users] Algorithm performance comparison with python

2016-03-08 Thread Jérémy Béjanin
On Tuesday, March 8, 2016 at 5:49:05 PM UTC-5, Stefan Karpinski wrote:

>
> You're writing the same matrix location each time so the results are the 
> same.
>

Once again you are correct, the problem was with the python code and I was 
confused because I had assumed it was right.

Back to the original problem:
I understand that this issue would not occur too often, considering that a 
serious algorithm would always use a sparse matrix (that's my next step!). 
Are there actual use cases where it might be important to use the same 
memory trick as Numpy?


Re: [julia-users] Algorithm performance comparison with python

2016-03-08 Thread Jérémy Béjanin

>
> We currently allocate an uninitialized array and then call the bzero 
> function to fill it with zeros – this takes a long time for large arrays. 
> We should instead mmap /dev/zero since that allows the kernel to lazily 
> allocate pages. I suspect this is what NumPy is doing, which in this 
> particular benchmark is an increasingly large benefit since huge chunks of 
> this matrix remain completely zeroed out and can all share a single actual 
> memory page. There's a very old (closed) issue about this that I've 
> reopened: #130 .
>

Thanks for looking into this.
 

> Also: I don't think this matters much but in the Julia version, the line
>
> ham[bra+1,bra+1] = diag_term
>
>
> is outside the `for s in 0:N-2` loop whereas in the Python version, the 
> corresponding line is inside the loop. That favors the Julia version, 
> however.
>

Yep, that was my mistake, thanks (strange I thought I had compared the 
outputs and they were the same, oh well...).


Re: [julia-users] Algorithm performance comparison with python

2016-03-08 Thread Jérémy Béjanin
I have. I don't see anything in particular that would help me here...
>
>

[julia-users] Algorithm performance comparison with python

2016-03-08 Thread Jérémy Béjanin
Hello,

By curiosity, I have translated a simple algorithm from python to julia in 
order to compare their performance (gist here 
). I am still a 
relative newcomer and so am not sure why I am seeing worse performance from 
julia.

The code was not optimized for python (or julia for that matters).

The code generates a Hamiltonian for a 1D transverse field ising model, and 
I am varying the number of sites which makes the size of the matrix grow 
exponentially. Here are the timings I get:



As can be seen, Julia is faster for lower number of sites, but that changes 
at N=13. Also interesting is the CPU and memory use during those runs. It 
seems that Julia uses memory much more aggressively:




[julia-users] Is build_executable.jl working?

2015-11-25 Thread Jérémy Béjanin
Hello,

I was wondering if the build_executable function was working. I tried using 
it on Windows and got the following error (I included some lines above the 
error for context):

...
> require.jl
> docs/helpdb.jl
> docs/basedocs.jl
> C:\Users\Jeremy\AppData\Local\Julia-0.4.1\share\julia\base\precompile.jl
> INFO: Linking sys.dll
> INFO: System image successfully built at 
> C:\Users\Jeremy\AppData\Local\Julia-0.4.1\bin\libtestitnow.dll
> WARNING: Building sys.dll on Windows against LLVM < 3.5.0 can cause 
> incorrect backtraces! Delete generated sys.dll to av
> oid these problems
> INFO: To run Julia with this image loaded, run: julia -J 
> C:\Users\Jeremy\AppData\Local\Julia-0.4.1\bin\libtestitnow.dll
>
> running: 
> C:\Users\Jeremy\.julia\v0.4\WinRPM\deps\usr\x86_64-w64-mingw32\sys-root\mingw\bin\gcc.exe
>  
> -g -Wl,--no-as-needed
> `-D_WIN32_WINNT=0x0502` 
> -IC:\Users\Jeremy\AppData\Local\Julia-0.4.1\include\julia 
> -IC:\Users\Jeremy\AppData\Local\src -
> IC:\Users\Jeremy\AppData\Local\src/support 
> -IC:\Users\Jeremy\AppData\Local\usr/include 
> -IC:\Users\Jeremy\.julia\v0.4\Win
> RPM\deps\usr\x86_64-w64-mingw32\sys-root\mingw\include 
> C:\Users\Jeremy\AppData\Local\Temp\jul12B2.tmp\start_func.c -o C:
> \Users\Jeremy\AppData\Local\Julia-0.4.1\bin\testitnow.exe 
> -Wl,-rpath,C:\Users\Jeremy\AppData\Local\Julia-0.4.1\bin -LC:\
> Users\Jeremy\AppData\Local\Julia-0.4.1\bin -ljulia -ltestitnow
> C:\Users\Jeremy\AppData\Local\Temp\jul12B2.tmp\start_func.c: In function 
> 'main':
> C:\Users\Jeremy\AppData\Local\Temp\jul12B2.tmp\start_func.c:19:5: error: 
> too few arguments to function 'jl_atexit_hook'
> jl_atexit_hook();
> ^
> In file included from 
> C:\Users\Jeremy\AppData\Local\Temp\jul12B2.tmp\start_func.c:1:0:
> C:\Users\Jeremy\AppData\Local\Julia-0.4.1\include\julia/julia.h:1188:16: 
> note: declared here
> DLLEXPORT void jl_atexit_hook(int status);
> ^
> ERROR: failed process: 
> Process(setenv(`'C:\Users\Jeremy\.julia\v0.4\WinRPM\deps\usr\x86_64-w64-mingw32\sys-root\mingw\bi
> n\gcc.exe' -g -Wl,--no-as-needed -D_WIN32_WINNT=0x0502 
> '-IC:\Users\Jeremy\AppData\Local\Julia-0.4.1\include\julia' '-IC:
> \Users\Jeremy\AppData\Local\src' 
> '-IC:\Users\Jeremy\AppData\Local\src/support' 
> '-IC:\Users\Jeremy\AppData\Local\usr/incl
> ude' 
> '-IC:\Users\Jeremy\.julia\v0.4\WinRPM\deps\usr\x86_64-w64-mingw32\sys-root\mingw\include'
>  
> 'C:\Users\Jeremy\AppData\
> Local\Temp\jul12B2.tmp\start_func.c' -o 
> 'C:\Users\Jeremy\AppData\Local\Julia-0.4.1\bin\testitnow.exe' 
> '-Wl,-rpath,C:\Use
> rs\Jeremy\AppData\Local\Julia-0.4.1\bin' 
> '-LC:\Users\Jeremy\AppData\Local\Julia-0.4.1\bin' -ljulia 
> -ltestitnow`,Union{AS
> CIIString,UTF8String}["=C:=C:\\Users\\Jeremy\\AppData\\Local\\Julia-0.4.1","ADMSXMLBINDIR=C:\\Program
>  
> Files (x86)\\Qucs\
>
> \bin","ADS_LICENSE_FILE=6609@localhost","ALLUSERSPROFILE=C:\\ProgramData","APPDATA=C:\\Users\\Jeremy\\AppData\\Roaming",
> "ASCOBINDIR=C:\\Program Files 
> (x86)\\Qucs\\bin","CommonProgramFiles=C:\\Program Files\\Common 
> Files","CommonProgramFiles
> (x86)=C:\\Program Files (x86)\\Common 
> Files","CommonProgramW6432=C:\\Program Files\\Common 
> Files","COMPUTERNAME=JBEJANIN
>
> ","ComSpec=C:\\Windows\\system32\\cmd.exe","FLEXLM_TIMEOUT=200","FP_NO_HOST_CHECK=NO","HOME=C:\\Users\\Jeremy","HOME
>
> DRIVE=C:","HOMEPATH=\\Users\\Jeremy","LM_LICENSE_FILE=6065@lka-ic-01;6065@localhost","LOCALAPPDATA=C:\\Users\\Jeremy\\Ap
> pData\\Local","LOGONSERVER=JBEJANIN","MATLAB_HOME=C:\\Program 
> Files\\MATLAB\\R2015a","MOZ_PLUGIN_PATH=C:\\PROGRAM FI
> LES (X86)\\FOXIT SOFTWARE\\FOXIT 
> READER\\plugins\\","NIDAQmxSwitchDir=C:\\Program Files (x86)\\National 
> Instruments\\NI-
> DAQ\\Switch\\","NUMBER_OF_PROCESSORS=4","OCTAVEDIR=C:\\Program Files 
> (x86)\\Qucs\\share\\qucs\\octave","OPENBLAS_MAIN_FR
>
> EE=1","OS=Windows_NT","Path=C:\\Users\\Jeremy\\AppData\\Local\\Julia-0.4.1\\bin;C:\\Users\\Jeremy\\AppData\\Local\\Julia
> -0.4.1\\bin\\..\\Git\\bin;C:\\Users\\Jeremy\\AppData\\Local\\Julia-0.4.1\\bin\\..\\Git\\usr\\bin;C:\\Program
>  
> Files (x86)
> \\SumatraPDF;C:\\Program Files (x86)\\NVIDIA 
> Corporation\\PhysX\\Common;C:\\ProgramData\\Oracle\\Java\\javapath;C:\\Prog
> ram Files 
> (x86)\\Qucs\\bin;C:\\Windows\\system32;C:\\Windows;C:\\Windows\\System32\\Wbem;C:\\Windows\\System32\\WindowsP
> owerShell\\v1.0\\;D:\\Programs\\MiKTeX 
> 2.9\\miktex\\bin\\x64\\;C:\\Users\\Jeremy\\AppData\\Local\\Julia-0.4.1\\bin;C:\\P
> rogram Files\\MATLAB\\R2015a\\runtime\\win64;C:\\Program 
> Files\\MATLAB\\R2015a\\bin;D:\\Programs\\Git\\cmd;D:\\Programs\
> \Git\\usr\\bin;D:\\Programs\\Vim\\vim74;C:\\Program Files 
> (x86)\\FastFieldSolvers\\FastHenry2\\Utilities;C:\\Program Fil
> es 
> (x86)\\FastFieldSolvers\\FastCap2\\Utilities;C:\\Users\\Jeremy\\AppData\\Local\\atom\\bin;;C:\\Users\\Jeremy\\.julia\
>
> \v0.4\\WinRPM\\deps\\usr\\x86_64-w64-mingw32\\sys-root\\mingw\\bin;C:\\Users\\Jeremy\\.julia\\v0.4\\WinRPM\\deps\\usr\\x
>
> 86_64-w64-mingw32\\sys-root\\mingw\\bin","PATHEXT

Re: [julia-users] Re: Reduction of multiple variables in parallel for-loop

2015-05-25 Thread Jérémy Béjanin
Thanks for the answers.

Would you have a better suggestion on how to do what I want?

As I explained, I am trying to fill two arrays where each array element is 
a function of the loop index. I realize that at the end of each loop, 
because the memory is not shared, each array will only have one element 
filled, and the rest will be zeroes. Adding all of these mostly zeroes 
arrays together to get the full one does not seem very efficient (or 
elegant).

Even if I had only one array to fill, the same efficiency issue exists.

I suppose a better way would be for each loop to return only the element, 
but I would also need to return each elements' array coordinate so that I 
could build the matrix properly at the end...

On Monday, May 25, 2015 at 2:57:38 PM UTC-4, Mohammed El-Beltagy wrote:
>
> Indeed. Also on second thought, I think that reduction on complex types is 
> not a very good idea. I suspect that with memory allocation at each step 
> and data copied across processes, it will make for a pretty slow parallel 
> code. 
>
> On Mon, May 25, 2015 at 8:15 PM, Toivo Henningsson  > wrote:
>
>> Though you shouldn't really redefine + for types that you haven't created 
>> yourself. But you could make your own container type where you define it, 
>> or maybe it's enough to put the two return values in an array. 
>
>
>

[julia-users] Reduction of multiple variables in parallel for-loop

2015-05-24 Thread Jérémy Béjanin
Is it possible to perform reduction on more than one variable in the case 
where the loop is computing more than one quantity?
I tried doing something like this:

var1=zeros(Float64,10,4)
var2=zeros(Complex128,10,1)
@parallel (+) for m=1:10
var1[m,:] = [m,2*m,3*m,4*m]
var2[m] = m
(var1, var2)
end

but the addition method cannot add that type...

The goal here is to fill two arrays where each loop computes a specific 
element. I figured that simply adding each end of loop array (which would 
not overlap since all non-modified entries are zeroes) would be more simple 
than using a shared array.

Thanks,
Jeremy


[julia-users] Pre-allocation and for-loops

2015-05-21 Thread Jérémy Béjanin
In the performance tips section of the manual, it is recommended to, if 
convenient/possible, pre-allocate the variables used in for-loops in order 
to avoid unecessary allocation and garbage collection.

The example given uses a function xinc(x) which returns a array . Running 
this function in a for loop results in a lot of memory allocation and takes 
significantly more time than the updating version of that same function xinc
!. This is clear to me.

However, if I change the line in xinc(x) to return a float: return [x, x+1, 
x+2] to return x, I notice that the memory use is only 96 bytes, which I 
assume means that no allocation is happening other than the initial one. 
There is no advantage in that case to using the xinc! function. I see that 
the manual suggests that the improvement will only occur with "complex" 
types. Why is that the case?

My other question concerns the necessity of updating functions: why is not 
possible to pre-allocate the variable and use the "normal" function, 
returning its output in the preallocated variable? The space in memory is 
already there, why is it not possible to use it? I believe that this is 
what MATLAB does.

Thanks for your help,
Jeremy